Line data Source code
1 : module radae
2 : !------------------------------------------------------------------------------
3 : !
4 : ! Description:
5 : !
6 : ! Data and subroutines to calculate absorptivities and emissivity needed
7 : ! for the LW radiation calculation.
8 : !
9 : ! Public interfaces are:
10 : !
11 : ! radae_init ------------ Initialization
12 : ! initialize_radbuffer -- Initialize the 3D abs/emis arrays.
13 : ! radabs ---------------- Compute absorptivities.
14 : ! radems ---------------- Compute emissivity.
15 : ! radtpl ---------------- Compute Temperatures and path lengths.
16 : ! radoz2 ---------------- Compute ozone path lengths.
17 : ! trcpth ---------------- Compute ghg path lengths.
18 : !
19 : ! Author: B. Collins
20 : !
21 : !------------------------------------------------------------------------------
22 : use shr_kind_mod, only: r8=>shr_kind_r8
23 : use spmd_utils, only: masterproc
24 : use ppgrid, only: pcols, pverp, begchunk, endchunk, pver
25 : use infnan, only: posinf, assignment(=)
26 : use pmgrid, only: plev, plevp
27 : use radconstants, only: nlwbands, idx_LW_0650_0800, idx_LW_0500_0650, &
28 : idx_LW_1000_1200, idx_LW_0800_1000, idx_LW_1200_2000
29 : use cam_abortutils, only: endrun
30 : use cam_logfile, only: iulog
31 : use wv_saturation, only: qsat_water
32 :
33 : implicit none
34 : private
35 : save
36 :
37 : public :: radabs, radems, radtpl, radae_init, radoz2, trcpth, initialize_radbuffer
38 :
39 : integer, public, parameter :: nbands = 2 ! Number of spectral bands
40 :
41 : ! Following data needed for restarts and in radclwmx
42 : real(r8), public, allocatable, target :: abstot_3d(:,:,:,:) ! Non-adjacent layer absorptivites
43 : real(r8), public, allocatable, target :: absnxt_3d(:,:,:,:) ! Nearest layer absorptivities
44 : real(r8), public, allocatable, target :: emstot_3d(:,:,:) ! Total emissivity
45 : integer, public :: ntoplw ! top level to solve for longwave cooling
46 :
47 : real(r8) :: p0 ! Standard pressure (dynes/cm**2)
48 : real(r8) :: amd ! Molecular weight of dry air (g/mol)
49 : real(r8) :: amco2 ! Molecular weight of co2 (g/mol)
50 : real(r8) :: mwo3 ! Molecular weight of O3 (g/mol)
51 :
52 : real(r8) :: gravit ! acceleration due to gravity (m/s**2)
53 : real(r8) :: gravit_cgs ! acceleration due to gravity (cm/s**2)
54 : real(r8) :: rga ! 1./gravit_cgs
55 : real(r8) :: epsilo ! Ratio of mol. wght of H2O to dry air
56 : real(r8) :: omeps ! 1._r8 - epsilo
57 : real(r8) :: sslp ! Standard sea-level pressure (dynes/cm**2)
58 : real(r8) :: stebol_cgs ! Stefan-Boltzmann's constant (CGS)
59 : real(r8) :: rgsslp ! 0.5/(gravit_cgs*sslp)
60 : real(r8) :: dpfo3 ! Voigt correction factor for O3
61 : real(r8) :: dpfco2 ! Voigt correction factor for CO2
62 :
63 : integer, parameter :: n_u = 25 ! Number of U in abs/emis tables
64 : integer, parameter :: n_p = 10 ! Number of P in abs/emis tables
65 : integer, parameter :: n_tp = 10 ! Number of T_p in abs/emis tables
66 : integer, parameter :: n_te = 21 ! Number of T_e in abs/emis tables
67 : integer, parameter :: n_rh = 7 ! Number of RH in abs/emis tables
68 :
69 : real(r8):: ah2onw(n_p, n_tp, n_u, n_te, n_rh) ! absorptivity (non-window)
70 : real(r8):: eh2onw(n_p, n_tp, n_u, n_te, n_rh) ! emissivity (non-window)
71 : real(r8):: ah2ow(n_p, n_tp, n_u, n_te, n_rh) ! absorptivity (window, for adjacent layers)
72 : real(r8):: cn_ah2ow(n_p, n_tp, n_u, n_te, n_rh) ! continuum transmission for absorptivity (window)
73 : real(r8):: cn_eh2ow(n_p, n_tp, n_u, n_te, n_rh) ! continuum transmission for emissivity (window)
74 : real(r8):: ln_ah2ow(n_p, n_tp, n_u, n_te, n_rh) ! line-only transmission for absorptivity (window)
75 : real(r8):: ln_eh2ow(n_p, n_tp, n_u, n_te, n_rh) ! line-only transmission for emissivity (window)
76 : !
77 : ! Constant coefficients for water vapor overlap with trace gases.
78 : ! Reference: Ramanathan, V. and P.Downey, 1986: A Nonisothermal
79 : ! Emissivity and Absorptivity Formulation for Water Vapor
80 : ! Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
81 : !
82 : real(r8):: coefh(2,4) = reshape( &
83 : (/ (/5.46557e+01_r8,-7.30387e-02_r8/), &
84 : (/1.09311e+02_r8,-1.46077e-01_r8/), &
85 : (/5.11479e+01_r8,-6.82615e-02_r8/), &
86 : (/1.02296e+02_r8,-1.36523e-01_r8/) /), (/2,4/) )
87 : !
88 : real(r8):: coefj(3,2) = reshape( &
89 : (/ (/2.82096e-02_r8,2.47836e-04_r8,1.16904e-06_r8/), &
90 : (/9.27379e-02_r8,8.04454e-04_r8,6.88844e-06_r8/) /), (/3,2/) )
91 : !
92 : real(r8):: coefk(3,2) = reshape( &
93 : (/ (/2.48852e-01_r8,2.09667e-03_r8,2.60377e-06_r8/) , &
94 : (/1.03594e+00_r8,6.58620e-03_r8,4.04456e-06_r8/) /), (/3,2/) )
95 : real(r8):: c16,c17,c26,c27,c28,c29,c30,c31
96 : !
97 : ! Farwing correction constants for narrow-band emissivity model,
98 : ! introduced to account for the deficiencies in narrow-band model
99 : ! used to derive the emissivity; tuned with Arkings line-by-line
100 : ! calculations. Just used for water vapor overlap with trace gases.
101 : !
102 : real(r8):: fwcoef ! Farwing correction constant
103 : real(r8):: fwc1,fwc2 ! Farwing correction constants
104 : real(r8):: fc1 ! Farwing correction constant
105 : !
106 : ! Collins/Hackney/Edwards (C/H/E) & Collins/Lee-Taylor/Edwards (C/LT/E)
107 : ! H2O parameterization
108 : !
109 : ! Notation:
110 : ! U = integral (P/P_0 dW) eq. 15 in Ramanathan/Downey 1986
111 : ! P = atmospheric pressure
112 : ! P_0 = reference atmospheric pressure
113 : ! W = precipitable water path
114 : ! T_e = emission temperature
115 : ! T_p = path temperature
116 : ! RH = path relative humidity
117 : !
118 : ! absorptivity/emissivity in window are fit using an expression:
119 : !
120 : ! a/e = f_a/e * {1.0 - ln_a/e * cn_a/e}
121 : !
122 : ! absorptivity/emissivity in non-window are fit using:
123 : !
124 : ! a/e = f_a/e * a/e_norm
125 : !
126 : ! where
127 : ! a/e = absorptivity/emissivity
128 : ! a/e_norm = absorptivity/emissivity normalized to 1
129 : ! f_a/e = value of a/e as U->infinity = f(T_e) only
130 : ! cn_a/e = continuum transmission
131 : ! ln_a/e = line transmission
132 : !
133 : ! spectral interval:
134 : ! 1 = 0-800 cm^-1 and 1200-2200 cm^-1 (rotation and rotation-vibration)
135 : ! 2 = 800-1200 cm^-1 (window)
136 : !
137 : real(r8), parameter:: min_tp_h2o = 160.0_r8 ! min T_p for pre-calculated abs/emis
138 : real(r8), parameter:: max_tp_h2o = 349.999999_r8 ! max T_p for pre-calculated abs/emis
139 : integer, parameter :: ntemp = 192 ! Number of temperatures in H2O sat. table for Tp
140 : integer, parameter :: o_fa = 6 ! Degree+1 of poly of T_e for absorptivity as U->inf.
141 : integer, parameter :: o_fe = 6 ! Degree+1 of poly of T_e for emissivity as U->inf.
142 : !-----------------------------------------------------------------------------
143 : ! Data for f in C/H/E fit -- value of A and E as U->infinity
144 : ! New C/LT/E fit (Hitran 2K, CKD 2.4) -- no change
145 : ! These values are determined by integrals of Planck functions or
146 : ! derivatives of Planck functions only.
147 : !-----------------------------------------------------------------------------
148 : !
149 : ! fa/fe coefficients for 2 bands (0-800 & 1200-2200, 800-1200 cm^-1)
150 : !
151 : ! Coefficients of polynomial for f_a in T_e
152 : !
153 : real(r8), parameter:: fat(o_fa,nbands) = reshape( (/ &
154 : (/-1.06665373E-01_r8, 2.90617375E-02_r8, -2.70642049E-04_r8, & ! 0-800&1200-2200 cm^-1
155 : 1.07595511E-06_r8, -1.97419681E-09_r8, 1.37763374E-12_r8/), & ! 0-800&1200-2200 cm^-1
156 : (/ 1.10666537E+00_r8, -2.90617375E-02_r8, 2.70642049E-04_r8, & ! 800-1200 cm^-1
157 : -1.07595511E-06_r8, 1.97419681E-09_r8, -1.37763374E-12_r8/) /) & ! 800-1200 cm^-1
158 : , (/o_fa,nbands/) )
159 : !
160 : ! Coefficients of polynomial for f_e in T_e
161 : !
162 : real(r8), parameter:: fet(o_fe,nbands) = reshape( (/ &
163 : (/3.46148163E-01_r8, 1.51240299E-02_r8, -1.21846479E-04_r8, & ! 0-800&1200-2200 cm^-1
164 : 4.04970123E-07_r8, -6.15368936E-10_r8, 3.52415071E-13_r8/), & ! 0-800&1200-2200 cm^-1
165 : (/6.53851837E-01_r8, -1.51240299E-02_r8, 1.21846479E-04_r8, & ! 800-1200 cm^-1
166 : -4.04970123E-07_r8, 6.15368936E-10_r8, -3.52415071E-13_r8/) /) & ! 800-1200 cm^-1
167 : , (/o_fa,nbands/) )
168 : !
169 : ! Note: max values should be slightly underestimated to avoid index bound violations
170 : !
171 : real(r8), parameter:: min_lp_h2o = -3.0_r8 ! min log_10(P) for pre-calculated abs/emis
172 : real(r8), parameter:: min_p_h2o = 1.0e-3_r8 ! min log_10(P) for pre-calculated abs/emis
173 : real(r8), parameter:: max_lp_h2o = -0.0000001_r8 ! max log_10(P) for pre-calculated abs/emis
174 : real(r8), parameter:: dlp_h2o = 0.3333333333333_r8 ! difference in adjacent elements of lp_h2o
175 :
176 : real(r8), parameter:: dtp_h2o = 21.111111111111_r8 ! difference in adjacent elements of tp_h2o
177 :
178 : real(r8), parameter:: min_rh_h2o = 0.0_r8 ! min RH for pre-calculated abs/emis
179 : real(r8), parameter:: max_rh_h2o = 1.19999999_r8 ! max RH for pre-calculated abs/emis
180 : real(r8), parameter:: drh_h2o = 0.2_r8 ! difference in adjacent elements of RH
181 :
182 : real(r8), parameter:: min_te_h2o = -120.0_r8 ! min T_e-T_p for pre-calculated abs/emis
183 : real(r8), parameter:: max_te_h2o = 79.999999_r8 ! max T_e-T_p for pre-calculated abs/emis
184 : real(r8), parameter:: dte_h2o = 10.0_r8 ! difference in adjacent elements of te_h2o
185 :
186 : real(r8), parameter:: min_lu_h2o = -8.0_r8 ! min log_10(U) for pre-calculated abs/emis
187 : real(r8), parameter:: min_u_h2o = 1.0e-8_r8 ! min pressure-weighted path-length
188 : real(r8), parameter:: max_lu_h2o = 3.9999999_r8 ! max log_10(U) for pre-calculated abs/emis
189 : real(r8), parameter:: dlu_h2o = 0.5_r8 ! difference in adjacent elements of lu_h2o
190 :
191 : real(r8), parameter:: g1(6)=(/0.0468556_r8,0.0397454_r8,0.0407664_r8,0.0304380_r8,0.0540398_r8,0.0321962_r8/)
192 : real(r8), parameter :: g2(6)=(/14.4832_r8,4.30242_r8,5.23523_r8,3.25342_r8,0.698935_r8,16.5599_r8/)
193 : real(r8), parameter :: g3(6)=(/26.1898_r8,18.4476_r8,15.3633_r8,12.1927_r8,9.14992_r8,8.07092_r8/)
194 : real(r8), parameter :: g4(6)=(/0.0261782_r8,0.0369516_r8,0.0307266_r8,0.0243854_r8,0.0182932_r8,0.0161418_r8/)
195 : real(r8), parameter :: ab(6)=(/3.0857e-2_r8,2.3524e-2_r8,1.7310e-2_r8,2.6661e-2_r8,2.8074e-2_r8,2.2915e-2_r8/)
196 : real(r8), parameter :: bb(6)=(/-1.3512e-4_r8,-6.8320e-5_r8,-3.2609e-5_r8,-1.0228e-5_r8,-9.5743e-5_r8,-1.0304e-4_r8/)
197 : real(r8), parameter :: abp(6)=(/2.9129e-2_r8,2.4101e-2_r8,1.9821e-2_r8,2.6904e-2_r8,2.9458e-2_r8,1.9892e-2_r8/)
198 : real(r8), parameter :: bbp(6)=(/-1.3139e-4_r8,-5.5688e-5_r8,-4.6380e-5_r8,-8.0362e-5_r8,-1.0115e-4_r8,-8.8061e-5_r8/)
199 :
200 :
201 : ! Public Interfaces
202 : !====================================================================================
203 : CONTAINS
204 : !====================================================================================
205 :
206 68112 : subroutine radabs(lchnk ,ncol , &
207 : pbr ,pnm ,co2em ,co2eml ,tplnka , &
208 : s2c ,tcg ,w ,h2otr ,plco2 , &
209 : plh2o ,co2t ,tint ,tlayr ,plol , &
210 : plos ,pmln ,piln ,ucfc11 ,ucfc12 , &
211 : un2o0 ,un2o1 ,uch4 ,uco211 ,uco212 , &
212 : uco213 ,uco221 ,uco222 ,uco223 ,uptype , &
213 : bn2o0 ,bn2o1 ,bch4 ,abplnk1 ,abplnk2 , &
214 68112 : abstot ,absnxt ,plh2ob ,wb , &
215 : odap_aer ,aer_trn_ttl, co2mmr)
216 : !-----------------------------------------------------------------------
217 : !
218 : ! Purpose:
219 : ! Compute absorptivities for h2o, co2, o3, ch4, n2o, cfc11 and cfc12
220 : !
221 : ! Method:
222 : ! h2o .... Uses nonisothermal emissivity method for water vapor from
223 : ! Ramanathan, V. and P.Downey, 1986: A Nonisothermal
224 : ! Emissivity and Absorptivity Formulation for Water Vapor
225 : ! Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
226 : !
227 : ! Implementation updated by Collins, Hackney, and Edwards (2001)
228 : ! using line-by-line calculations based upon Hitran 1996 and
229 : ! CKD 2.1 for absorptivity and emissivity
230 : !
231 : ! Implementation updated by Collins, Lee-Taylor, and Edwards (2003)
232 : ! using line-by-line calculations based upon Hitran 2000 and
233 : ! CKD 2.4 for absorptivity and emissivity
234 : !
235 : ! co2 .... Uses absorptance parameterization of the 15 micro-meter
236 : ! (500 - 800 cm-1) band system of Carbon Dioxide, from
237 : ! Kiehl, J.T. and B.P.Briegleb, 1991: A New Parameterization
238 : ! of the Absorptance Due to the 15 micro-meter Band System
239 : ! of Carbon Dioxide Jouranl of Geophysical Research,
240 : ! vol. 96., D5, pp 9013-9019.
241 : ! Parameterizations for the 9.4 and 10.4 mircon bands of CO2
242 : ! are also included.
243 : !
244 : ! o3 .... Uses absorptance parameterization of the 9.6 micro-meter
245 : ! band system of ozone, from Ramanathan, V. and R.Dickinson,
246 : ! 1979: The Role of stratospheric ozone in the zonal and
247 : ! seasonal radiative energy balance of the earth-troposphere
248 : ! system. Journal of the Atmospheric Sciences, Vol. 36,
249 : ! pp 1084-1104
250 : !
251 : ! ch4 .... Uses a broad band model for the 7.7 micron band of methane.
252 : !
253 : ! n20 .... Uses a broad band model for the 7.8, 8.6 and 17.0 micron
254 : ! bands of nitrous oxide
255 : !
256 : ! cfc11 ... Uses a quasi-linear model for the 9.2, 10.7, 11.8 and 12.5
257 : ! micron bands of CFC11
258 : !
259 : ! cfc12 ... Uses a quasi-linear model for the 8.6, 9.1, 10.8 and 11.2
260 : ! micron bands of CFC12
261 : !
262 : !
263 : ! Computes individual absorptivities for non-adjacent layers, accounting
264 : ! for band overlap, and sums to obtain the total; then, computes the
265 : ! nearest layer contribution.
266 : !
267 : ! Author: W. Collins (H2O absorptivity) and J. Kiehl
268 : !------------------------------Arguments--------------------------------
269 : !
270 : ! Input arguments
271 : !
272 : integer, intent(in) :: lchnk ! chunk identifier
273 : integer, intent(in) :: ncol ! number of atmospheric columns
274 :
275 : real(r8), intent(in) :: pbr(pcols,pver) ! Prssr at mid-levels (dynes/cm2)
276 : real(r8), intent(in) :: pnm(pcols,pverp) ! Prssr at interfaces (dynes/cm2)
277 : real(r8), intent(in) :: co2em(pcols,pverp) ! Co2 emissivity function
278 : real(r8), intent(in) :: co2eml(pcols,pver) ! Co2 emissivity function
279 : real(r8), intent(in) :: tplnka(pcols,pverp) ! Planck fnctn level temperature
280 : real(r8), intent(in) :: s2c(pcols,pverp) ! H2o continuum path length
281 : real(r8), intent(in) :: tcg(pcols,pverp) ! H2o-mass-wgted temp. (Curtis-Godson approx.)
282 : real(r8), intent(in) :: w(pcols,pverp) ! H2o prs wghted path
283 : real(r8), intent(in) :: h2otr(pcols,pverp) ! H2o trnsmssn fnct for o3 overlap
284 : real(r8), intent(in) :: plco2(pcols,pverp) ! Co2 prs wghted path length
285 : real(r8), intent(in) :: plh2o(pcols,pverp) ! H2o prs wfhted path length
286 : real(r8), intent(in) :: co2t(pcols,pverp) ! Tmp and prs wghted path length
287 : real(r8), intent(in) :: tint(pcols,pverp) ! Interface temperatures
288 : real(r8), intent(in) :: tlayr(pcols,pverp) ! K-1 level temperatures
289 : real(r8), intent(in) :: plol(pcols,pverp) ! Ozone prs wghted path length
290 : real(r8), intent(in) :: plos(pcols,pverp) ! Ozone path length
291 : real(r8), intent(in) :: pmln(pcols,pver) ! Ln(pmidm1)
292 : real(r8), intent(in) :: piln(pcols,pverp) ! Ln(pintm1)
293 : real(r8), intent(in) :: plh2ob(nbands,pcols,pverp) ! Pressure weighted h2o path with
294 : ! Hulst-Curtis-Godson temp. factor
295 : ! for H2O bands
296 : real(r8), intent(in) :: wb(nbands,pcols,pverp) ! H2o path length with
297 : ! Hulst-Curtis-Godson temp. factor
298 : ! for H2O bands
299 :
300 : ! [fraction] absorbtion optical depth, cumulative from top
301 : real(r8), intent(in) :: odap_aer(pcols,pver,nlwbands)
302 :
303 : ! [fraction] Total transmission between interfaces k1 and k2
304 : real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands)
305 :
306 : !
307 : ! Trace gas variables
308 : !
309 : real(r8), intent(in) :: co2mmr(pcols) ! co2 column mean mass mixing ratio
310 : real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
311 : real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
312 : real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
313 : real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
314 : real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
315 : real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
316 : real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
317 : real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
318 : real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
319 : real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
320 : real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
321 : real(r8), intent(in) :: uptype(pcols,pverp) ! continuum path length
322 : real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
323 : real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
324 : real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
325 : real(r8), intent(in) :: abplnk1(14,pcols,pverp) ! non-nearest layer Planck factor
326 : real(r8), intent(in) :: abplnk2(14,pcols,pverp) ! nearest layer factor
327 : !
328 : ! Output arguments
329 : !
330 : real(r8), intent(out) :: abstot(pcols,ntoplw:pverp,ntoplw:pverp) ! Total absorptivity
331 : real(r8), intent(out) :: absnxt(pcols,pver,4) ! Total nearest layer absorptivity
332 : !
333 : !---------------------------Local variables-----------------------------
334 : !
335 : integer i ! Longitude index
336 : integer k ! Level index
337 : integer k1 ! Level index
338 : integer k2 ! Level index
339 : integer kn ! Nearest level index
340 : integer wvl ! Wavelength index
341 :
342 : real(r8) abstrc(pcols) ! total trace gas absorptivity
343 : real(r8) bplnk(14,pcols,4) ! Planck functions for sub-divided layers
344 : real(r8) pnew(pcols) ! Effective pressure for H2O vapor linewidth
345 : real(r8) pnewb(nbands) ! Effective pressure for h2o linewidth w/
346 : ! Hulst-Curtis-Godson correction for
347 : ! each band
348 : real(r8) u(pcols) ! Pressure weighted H2O path length
349 : real(r8) ub(nbands) ! Pressure weighted H2O path length with
350 : ! Hulst-Curtis-Godson correction for
351 : ! each band
352 : real(r8) tbar(pcols,4) ! Mean layer temperature
353 : real(r8) emm(pcols,4) ! Mean co2 emissivity
354 : real(r8) o3emm(pcols,4) ! Mean o3 emissivity
355 : real(r8) o3bndi ! Ozone band parameter
356 : real(r8) temh2o(pcols,4) ! Mean layer temperature equivalent to tbar
357 : real(r8) k21 ! Exponential coefficient used to calculate
358 : ! ! rotation band transmissvty in the 650-800
359 : ! ! cm-1 region (tr1)
360 : real(r8) k22 ! Exponential coefficient used to calculate
361 : ! ! rotation band transmissvty in the 500-650
362 : ! ! cm-1 region (tr2)
363 : real(r8) uc1(pcols) ! H2o continuum pathlength in 500-800 cm-1
364 : real(r8) to3h2o(pcols) ! H2o trnsmsn for overlap with o3
365 : real(r8) pi ! For co2 absorptivity computation
366 : real(r8) sqti(pcols) ! Used to store sqrt of mean temperature
367 : real(r8) et ! Co2 hot band factor
368 : real(r8) et2 ! Co2 hot band factor squared
369 : real(r8) et4 ! Co2 hot band factor to fourth power
370 : real(r8) omet ! Co2 stimulated emission term
371 : real(r8) f1co2 ! Co2 central band factor
372 : real(r8) f2co2(pcols) ! Co2 weak band factor
373 : real(r8) f3co2(pcols) ! Co2 weak band factor
374 : real(r8) t1co2(pcols) ! Overlap factr weak bands on strong band
375 : real(r8) sqwp ! Sqrt of co2 pathlength
376 : real(r8) f1sqwp(pcols) ! Main co2 band factor
377 : real(r8) oneme ! Co2 stimulated emission term
378 : real(r8) alphat ! Part of the co2 stimulated emission term
379 : real(r8) co2vmr(pcols) ! CO2 column mean vmr
380 : real(r8) rmw ! ratio of molecular weights (air/co2)
381 : real(r8) wco2 ! Constants used to define co2 pathlength
382 : real(r8) posqt ! Effective pressure for co2 line width
383 : real(r8) u7(pcols) ! Co2 hot band path length
384 : real(r8) u8 ! Co2 hot band path length
385 : real(r8) u9 ! Co2 hot band path length
386 : real(r8) u13 ! Co2 hot band path length
387 : real(r8) rbeta7(pcols) ! Inverse of co2 hot band line width par
388 : real(r8) rbeta8 ! Inverse of co2 hot band line width par
389 : real(r8) rbeta9 ! Inverse of co2 hot band line width par
390 : real(r8) rbeta13 ! Inverse of co2 hot band line width par
391 : real(r8) tpatha ! For absorptivity computation
392 : real(r8) abso(pcols,4) ! Absorptivity for various gases/bands
393 : real(r8) dtx(pcols) ! Planck temperature minus 250 K
394 : real(r8) dty(pcols) ! Path temperature minus 250 K
395 : real(r8) term7(pcols,2) ! Kl_inf(i) in eq(r8) of table A3a of R&D
396 : real(r8) term8(pcols,2) ! Delta kl_inf(i) in eq(r8)
397 : real(r8) tr1 ! Eqn(6) in table A2 of R&D for 650-800
398 : real(r8) tr10(pcols) ! Eqn (6) times eq(4) in table A2
399 : ! ! of R&D for 500-650 cm-1 region
400 : real(r8) tr2 ! Eqn(6) in table A2 of R&D for 500-650
401 : real(r8) tr5 ! Eqn(4) in table A2 of R&D for 650-800
402 : real(r8) tr6 ! Eqn(4) in table A2 of R&D for 500-650
403 : real(r8) tr9(pcols) ! Equation (6) times eq(4) in table A2
404 : ! ! of R&D for 650-800 cm-1 region
405 : real(r8) sqrtu(pcols) ! Sqrt of pressure weighted h20 pathlength
406 : real(r8) fwk(pcols) ! Equation(33) in R&D far wing correction
407 : real(r8) fwku(pcols) ! GU term in eqs(1) and (6) in table A2
408 : real(r8) to3co2(pcols) ! P weighted temp in ozone band model
409 : real(r8) dpnm(pcols) ! Pressure difference between two levels
410 : real(r8) pnmsq(pcols,pverp) ! Pressure squared
411 : real(r8) dw(pcols) ! Amount of h2o between two levels
412 : real(r8) uinpl(pcols,4) ! Nearest layer subdivision factor
413 : real(r8) winpl(pcols,4) ! Nearest layer subdivision factor
414 : real(r8) zinpl(pcols,4) ! Nearest layer subdivision factor
415 : real(r8) pinpl(pcols,4) ! Nearest layer subdivision factor
416 : real(r8) dplh2o(pcols) ! Difference in press weighted h2o amount
417 : real(r8) r293 ! 1/293
418 : real(r8) r250 ! 1/250
419 : real(r8) r3205 ! Line width factor for o3 (see R&Di)
420 : real(r8) r300 ! 1/300
421 : real(r8) rsslp ! Reciprocal of sea level pressure
422 : real(r8) r2sslp ! 1/2 of rsslp
423 : real(r8) ds2c ! Y in eq(7) in table A2 of R&D
424 : real(r8) dplos ! Ozone pathlength eq(A2) in R&Di
425 : real(r8) dplol ! Presure weighted ozone pathlength
426 : real(r8) tlocal ! Local interface temperature
427 : real(r8) beta ! Ozone mean line parameter eq(A3) in R&Di
428 : ! (includes Voigt line correction factor)
429 : real(r8) rphat ! Effective pressure for ozone beta
430 : real(r8) tcrfac ! Ozone temperature factor table 1 R&Di
431 : real(r8) tmp1 ! Ozone band factor see eq(A1) in R&Di
432 : real(r8) u1 ! Effective ozone pathlength eq(A2) in R&Di
433 : real(r8) realnu ! 1/beta factor in ozone band model eq(A1)
434 : real(r8) tmp2 ! Ozone band factor see eq(A1) in R&Di
435 : real(r8) u2 ! Effective ozone pathlength eq(A2) in R&Di
436 : real(r8) rsqti ! Reciprocal of sqrt of path temperature
437 : real(r8) tpath ! Path temperature used in co2 band model
438 : real(r8) tmp3 ! Weak band factor see K&B
439 : real(r8) rdpnmsq ! Reciprocal of difference in press^2
440 : real(r8) rdpnm ! Reciprocal of difference in press
441 : real(r8) p1 ! Mean pressure factor
442 : real(r8) p2 ! Mean pressure factor
443 : real(r8) dtym10 ! T - 260 used in eq(9) and (10) table A3a
444 : real(r8) dplco2 ! Co2 path length
445 : real(r8) te ! A_0 T factor in ozone model table 1 of R&Di
446 : real(r8) denom ! Denominator in eq(r8) of table A3a of R&D
447 : real(r8) th2o(pcols) ! transmission due to H2O
448 : real(r8) tco2(pcols) ! transmission due to CO2
449 : real(r8) to3(pcols) ! transmission due to O3
450 : !
451 : ! Transmission terms for various spectral intervals:
452 : !
453 : real(r8) trab2(pcols) ! H2o 500 - 800 cm-1
454 : real(r8) absbnd ! Proportional to co2 band absorptance
455 : real(r8) dbvtit(pcols,pverp)! Intrfc drvtv plnck fnctn for o3
456 : real(r8) dbvtly(pcols,pver) ! Level drvtv plnck fnctn for o3
457 : !
458 : ! Variables for Collins/Hackney/Edwards (C/H/E) &
459 : ! Collins/Lee-Taylor/Edwards (C/LT/E) H2O parameterization
460 :
461 : !
462 : ! Notation:
463 : ! U = integral (P/P_0 dW) eq. 15 in Ramanathan/Downey 1986
464 : ! P = atmospheric pressure
465 : ! P_0 = reference atmospheric pressure
466 : ! W = precipitable water path
467 : ! T_e = emission temperature
468 : ! T_p = path temperature
469 : ! RH = path relative humidity
470 : !
471 : real(r8) fa ! asymptotic value of abs. as U->infinity
472 : real(r8) a_star ! normalized absorptivity for non-window
473 : real(r8) l_star ! interpolated line transmission
474 : real(r8) c_star ! interpolated continuum transmission
475 :
476 : real(r8) te1 ! emission temperature
477 : real(r8) te2 ! te^2
478 : real(r8) te3 ! te^3
479 : real(r8) te4 ! te^4
480 : real(r8) te5 ! te^5
481 :
482 : real(r8) log_u ! log base 10 of U
483 : real(r8) log_uc ! log base 10 of H2O continuum path
484 : real(r8) log_p ! log base 10 of P
485 : real(r8) t_p ! T_p
486 : real(r8) t_e ! T_e (offset by T_p)
487 :
488 : integer iu ! index for log10(U)
489 : integer iu1 ! iu + 1
490 : integer iuc ! index for log10(H2O continuum path)
491 : integer iuc1 ! iuc + 1
492 : integer ip ! index for log10(P)
493 : integer ip1 ! ip + 1
494 : integer itp ! index for T_p
495 : integer itp1 ! itp + 1
496 : integer ite ! index for T_e
497 : integer ite1 ! ite + 1
498 : integer irh ! index for RH
499 : integer irh1 ! irh + 1
500 :
501 : real(r8) dvar ! normalized variation in T_p/T_e/P/U
502 : real(r8) uvar ! U * diffusivity factor
503 : real(r8) uscl ! factor for lineary scaling as U->0
504 :
505 : real(r8) wu ! weight for U
506 : real(r8) wu1 ! 1 - wu
507 : real(r8) wuc ! weight for H2O continuum path
508 : real(r8) wuc1 ! 1 - wuc
509 : real(r8) wp ! weight for P
510 : real(r8) wp1 ! 1 - wp
511 : real(r8) wtp ! weight for T_p
512 : real(r8) wtp1 ! 1 - wtp
513 : real(r8) wte ! weight for T_e
514 : real(r8) wte1 ! 1 - wte
515 : real(r8) wrh ! weight for RH
516 : real(r8) wrh1 ! 1 - wrh
517 :
518 : real(r8) w_0_0_ ! weight for Tp/Te combination
519 : real(r8) w_0_1_ ! weight for Tp/Te combination
520 : real(r8) w_1_0_ ! weight for Tp/Te combination
521 : real(r8) w_1_1_ ! weight for Tp/Te combination
522 :
523 : real(r8) w_0_00 ! weight for Tp/Te/RH combination
524 : real(r8) w_0_01 ! weight for Tp/Te/RH combination
525 : real(r8) w_0_10 ! weight for Tp/Te/RH combination
526 : real(r8) w_0_11 ! weight for Tp/Te/RH combination
527 : real(r8) w_1_00 ! weight for Tp/Te/RH combination
528 : real(r8) w_1_01 ! weight for Tp/Te/RH combination
529 : real(r8) w_1_10 ! weight for Tp/Te/RH combination
530 : real(r8) w_1_11 ! weight for Tp/Te/RH combination
531 :
532 : real(r8) w00_00 ! weight for P/Tp/Te/RH combination
533 : real(r8) w00_01 ! weight for P/Tp/Te/RH combination
534 : real(r8) w00_10 ! weight for P/Tp/Te/RH combination
535 : real(r8) w00_11 ! weight for P/Tp/Te/RH combination
536 : real(r8) w01_00 ! weight for P/Tp/Te/RH combination
537 : real(r8) w01_01 ! weight for P/Tp/Te/RH combination
538 : real(r8) w01_10 ! weight for P/Tp/Te/RH combination
539 : real(r8) w01_11 ! weight for P/Tp/Te/RH combination
540 : real(r8) w10_00 ! weight for P/Tp/Te/RH combination
541 : real(r8) w10_01 ! weight for P/Tp/Te/RH combination
542 : real(r8) w10_10 ! weight for P/Tp/Te/RH combination
543 : real(r8) w10_11 ! weight for P/Tp/Te/RH combination
544 : real(r8) w11_00 ! weight for P/Tp/Te/RH combination
545 : real(r8) w11_01 ! weight for P/Tp/Te/RH combination
546 : real(r8) w11_10 ! weight for P/Tp/Te/RH combination
547 : real(r8) w11_11 ! weight for P/Tp/Te/RH combination
548 :
549 : integer ib ! spectral interval:
550 : ! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
551 : ! 2 = 800-1200 cm^-1
552 :
553 :
554 : real(r8) pch2o ! H2O continuum path
555 : real(r8) fch2o ! temp. factor for continuum
556 : real(r8) uch2o ! U corresponding to H2O cont. path (window)
557 :
558 : real(r8) fdif ! secant(zenith angle) for diffusivity approx.
559 :
560 : real(r8) sslp_mks ! Sea-level pressure in MKS units
561 : real(r8) esx ! saturation vapor pressure returned by qsat
562 : real(r8) qsx ! saturation mixing ratio returned by qsat
563 : real(r8) pnew_mks ! pnew in MKS units
564 : real(r8) q_path ! effective specific humidity along path
565 : real(r8) rh_path ! effective relative humidity along path
566 :
567 : integer bnd_idx ! LW band index
568 : real(r8) aer_pth_dlt ! [kg m-2] STRAER path between interface levels k1 and k2
569 : real(r8) aer_pth_ngh(pcols)
570 : ! [kg m-2] STRAER path between neighboring layers
571 : real(r8) odap_aer_ttl ! [fraction] Total path absorption optical depth
572 : real(r8) aer_trn_ngh(pcols,nlwbands)
573 : ! [fraction] Total transmission between
574 : ! nearest neighbor sub-levels
575 : !
576 : !--------------------------Statement function---------------------------
577 : !
578 : real(r8) dbvt,t ! Planck fnctn tmp derivative for o3
579 : !
580 : dbvt(t)=(-2.8911366682e-4_r8+(2.3771251896e-6_r8+1.1305188929e-10_r8*t)*t)/ &
581 : (1.0_r8+(-6.1364820707e-3_r8+1.5550319767e-5_r8*t)*t)
582 : !
583 : !
584 : !-----------------------------------------------------------------------
585 : !
586 : ! Initialize
587 : !
588 340560 : do k2=1,4
589 340560 : do k1=1,ntoplw-1
590 272448 : absnxt(:,k1,k2) = posinf ! set unused portions for lf95 restart write
591 : end do
592 : end do
593 :
594 1907136 : do k=ntoplw,pverp
595 1907136 : abstot(:,k,k) = posinf ! set unused portions for lf95 restart write
596 : end do
597 :
598 1839024 : do k=ntoplw,pver
599 29638224 : do i=1,ncol
600 27799200 : dbvtly(i,k) = dbvt(tlayr(i,k+1))
601 29570112 : dbvtit(i,k) = dbvt(tint(i,k))
602 : end do
603 : end do
604 68112 : rmw = amd/amco2
605 1137312 : do i=1,ncol
606 1069200 : dbvtit(i,pverp) = dbvt(tint(i,pverp))
607 1137312 : co2vmr(i) = co2mmr(i) * rmw
608 : end do
609 : !
610 68112 : r293 = 1._r8/293._r8
611 68112 : r250 = 1._r8/250._r8
612 68112 : r3205 = 1._r8/.3205_r8
613 68112 : r300 = 1._r8/300._r8
614 68112 : rsslp = 1._r8/sslp
615 68112 : r2sslp = 1._r8/(2._r8*sslp)
616 : !
617 : !Constants for computing U corresponding to H2O cont. path
618 : !
619 68112 : fdif = 1.66_r8
620 68112 : sslp_mks = sslp / 10.0_r8
621 : !
622 : ! Non-adjacent layer absorptivity:
623 : !
624 : ! abso(i,1) 0 - 800 cm-1 h2o rotation band
625 : ! abso(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
626 : ! abso(i,2) 800 - 1200 cm-1 h2o window
627 : !
628 : ! Separation between rotation and vibration-rotation dropped, so
629 : ! only 2 slots needed for H2O absorptivity
630 : !
631 : ! 500-800 cm^-1 H2o continuum/line overlap already included
632 : ! in abso(i,1). This used to be in abso(i,4)
633 : !
634 : ! abso(i,3) o3 9.6 micrometer band (nu3 and nu1 bands)
635 : ! abso(i,4) co2 15 micrometer band system
636 : !
637 1907136 : do k=ntoplw,pverp
638 30775536 : do i=1,ncol
639 28868400 : pnmsq(i,k) = pnm(i,k)**2
640 30707424 : dtx(i) = tplnka(i,k) - 250._r8
641 : end do
642 : end do
643 : !
644 : ! Non-nearest layer level loops
645 : !
646 1907136 : do k1=pverp,ntoplw,-1
647 51560784 : do k2=pverp,ntoplw,-1
648 49653648 : if (k1 == k2) cycle
649 798393024 : do i=1,ncol
650 750578400 : dplh2o(i) = plh2o(i,k1) - plh2o(i,k2)
651 750578400 : u(i) = abs(dplh2o(i))
652 750578400 : sqrtu(i) = sqrt(u(i))
653 750578400 : ds2c = abs(s2c(i,k1) - s2c(i,k2))
654 750578400 : dw(i) = abs(w(i,k1) - w(i,k2))
655 750578400 : uc1(i) = (ds2c + 1.7e-3_r8*u(i))*(1._r8 + 2._r8*ds2c)/(1._r8 + 15._r8*ds2c)
656 750578400 : pch2o = ds2c
657 750578400 : pnew(i) = u(i)/dw(i)
658 750578400 : pnew_mks = pnew(i) * sslp_mks
659 : !
660 : ! Changed effective path temperature to std. Curtis-Godson form
661 : !
662 750578400 : tpatha = abs(tcg(i,k1) - tcg(i,k2))/dw(i)
663 750578400 : t_p = min(max(tpatha, min_tp_h2o), max_tp_h2o)
664 :
665 750578400 : call qsat_water(t_p, pnew_mks, esx, qsx)
666 : !
667 : ! Compute effective RH along path
668 : !
669 750578400 : q_path = dw(i) / abs(pnm(i,k1) - pnm(i,k2)) / rga
670 : !
671 : ! Calculate effective u, pnew for each band using
672 : ! Hulst-Curtis-Godson approximation:
673 : ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
674 : ! 2nd edition, Oxford University Press, 1989.
675 : ! Effective H2O path (w)
676 : ! eq. 6.24, p. 228
677 : ! Effective H2O path pressure (pnew = u/w):
678 : ! eq. 6.29, p. 228
679 : !
680 750578400 : ub(1) = abs(plh2ob(1,i,k1) - plh2ob(1,i,k2)) / psi(t_p,1)
681 750578400 : ub(2) = abs(plh2ob(2,i,k1) - plh2ob(2,i,k2)) / psi(t_p,2)
682 :
683 750578400 : pnewb(1) = ub(1) / abs(wb(1,i,k1) - wb(1,i,k2)) * phi(t_p,1)
684 750578400 : pnewb(2) = ub(2) / abs(wb(2,i,k1) - wb(2,i,k2)) * phi(t_p,2)
685 :
686 750578400 : dtx(i) = tplnka(i,k2) - 250._r8
687 750578400 : dty(i) = tpatha - 250._r8
688 :
689 750578400 : fwk(i) = fwcoef + fwc1/(1._r8 + fwc2*u(i))
690 750578400 : fwku(i) = fwk(i)*u(i)
691 : !
692 : ! Define variables for C/H/E (now C/LT/E) fit
693 : !
694 : ! abso(i,1) 0 - 800 cm-1 h2o rotation band
695 : ! abso(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
696 : ! abso(i,2) 800 - 1200 cm-1 h2o window
697 : !
698 : ! Separation between rotation and vibration-rotation dropped, so
699 : ! only 2 slots needed for H2O absorptivity
700 : !
701 : ! Notation:
702 : ! U = integral (P/P_0 dW)
703 : ! P = atmospheric pressure
704 : ! P_0 = reference atmospheric pressure
705 : ! W = precipitable water path
706 : ! T_e = emission temperature
707 : ! T_p = path temperature
708 : ! RH = path relative humidity
709 : !
710 : !
711 : ! Terms for asymptotic value of emissivity
712 : !
713 750578400 : te1 = tplnka(i,k2)
714 750578400 : te2 = te1 * te1
715 750578400 : te3 = te2 * te1
716 750578400 : te4 = te3 * te1
717 750578400 : te5 = te4 * te1
718 :
719 : !
720 : ! Band-independent indices for lines and continuum tables
721 : !
722 750578400 : dvar = (t_p - min_tp_h2o) / dtp_h2o
723 750578400 : itp = min(max(int(aint(dvar,r8)) + 1, 1), n_tp - 1)
724 750578400 : itp1 = itp + 1
725 750578400 : wtp = dvar - floor(dvar)
726 750578400 : wtp1 = 1.0_r8 - wtp
727 :
728 750578400 : t_e = min(max(tplnka(i,k2)-t_p, min_te_h2o), max_te_h2o)
729 750578400 : dvar = (t_e - min_te_h2o) / dte_h2o
730 750578400 : ite = min(max(int(aint(dvar,r8)) + 1, 1), n_te - 1)
731 750578400 : ite1 = ite + 1
732 750578400 : wte = dvar - floor(dvar)
733 750578400 : wte1 = 1.0_r8 - wte
734 :
735 750578400 : rh_path = min(max(q_path / qsx, min_rh_h2o), max_rh_h2o)
736 750578400 : dvar = (rh_path - min_rh_h2o) / drh_h2o
737 750578400 : irh = min(max(int(aint(dvar,r8)) + 1, 1), n_rh - 1)
738 750578400 : irh1 = irh + 1
739 750578400 : wrh = dvar - floor(dvar)
740 750578400 : wrh1 = 1.0_r8 - wrh
741 :
742 750578400 : w_0_0_ = wtp * wte
743 750578400 : w_0_1_ = wtp * wte1
744 750578400 : w_1_0_ = wtp1 * wte
745 750578400 : w_1_1_ = wtp1 * wte1
746 :
747 750578400 : w_0_00 = w_0_0_ * wrh
748 750578400 : w_0_01 = w_0_0_ * wrh1
749 750578400 : w_0_10 = w_0_1_ * wrh
750 750578400 : w_0_11 = w_0_1_ * wrh1
751 750578400 : w_1_00 = w_1_0_ * wrh
752 750578400 : w_1_01 = w_1_0_ * wrh1
753 750578400 : w_1_10 = w_1_1_ * wrh
754 750578400 : w_1_11 = w_1_1_ * wrh1
755 :
756 : !
757 : ! H2O Continuum path for 0-800 and 1200-2200 cm^-1
758 : !
759 : ! Assume foreign continuum dominates total H2O continuum in these bands
760 : ! per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776
761 : ! Then the effective H2O path is just
762 : ! U_c = integral[ f(P) dW ]
763 : ! where
764 : ! W = water-vapor mass and
765 : ! f(P) = dependence of foreign continuum on pressure
766 : ! = P / sslp
767 : ! Then
768 : ! U_c = U (the same effective H2O path as for lines)
769 : !
770 : !
771 : ! Continuum terms for 800-1200 cm^-1
772 : !
773 : ! Assume self continuum dominates total H2O continuum for this band
774 : ! per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776
775 : ! Then the effective H2O self-continuum path is
776 : ! U_c = integral[ h(e,T) dW ] (*eq. 1*)
777 : ! where
778 : ! W = water-vapor mass and
779 : ! e = partial pressure of H2O along path
780 : ! T = temperature along path
781 : ! h(e,T) = dependence of foreign continuum on e,T
782 : ! = e / sslp * f(T)
783 : !
784 : ! Replacing
785 : ! e =~ q * P / epsilo
786 : ! q = mixing ratio of H2O
787 : ! epsilo = 0.622
788 : !
789 : ! and using the definition
790 : ! U = integral [ (P / sslp) dW ]
791 : ! = (P / sslp) W (homogeneous path)
792 : !
793 : ! the effective path length for the self continuum is
794 : ! U_c = (q / epsilo) f(T) U (*eq. 2*)
795 : !
796 : ! Once values of T, U, and q have been calculated for the inhomogeneous
797 : ! path, this sets U_c for the corresponding
798 : ! homogeneous atmosphere. However, this need not equal the
799 : ! value of U_c' defined by eq. 1 for the actual inhomogeneous atmosphere
800 : ! under consideration.
801 : !
802 : ! Solution: hold T and q constant, solve for U' that gives U_c' by
803 : ! inverting eq. (2):
804 : !
805 : ! U' = (U_c * epsilo) / (q * f(T))
806 : !
807 750578400 : fch2o = fh2oself(t_p)
808 750578400 : uch2o = (pch2o * epsilo) / (q_path * fch2o)
809 :
810 : !
811 : ! Band-dependent indices for non-window
812 : !
813 750578400 : ib = 1
814 :
815 750578400 : uvar = ub(ib) * fdif
816 750578400 : log_u = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
817 750578400 : dvar = (log_u - min_lu_h2o) / dlu_h2o
818 750578400 : iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
819 750578400 : iu1 = iu + 1
820 750578400 : wu = dvar - floor(dvar)
821 750578400 : wu1 = 1.0_r8 - wu
822 :
823 750578400 : log_p = min(log10(max(pnewb(ib), min_p_h2o)), max_lp_h2o)
824 750578400 : dvar = (log_p - min_lp_h2o) / dlp_h2o
825 750578400 : ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
826 750578400 : ip1 = ip + 1
827 750578400 : wp = dvar - floor(dvar)
828 750578400 : wp1 = 1.0_r8 - wp
829 :
830 750578400 : w00_00 = wp * w_0_00
831 750578400 : w00_01 = wp * w_0_01
832 750578400 : w00_10 = wp * w_0_10
833 750578400 : w00_11 = wp * w_0_11
834 750578400 : w01_00 = wp * w_1_00
835 750578400 : w01_01 = wp * w_1_01
836 750578400 : w01_10 = wp * w_1_10
837 750578400 : w01_11 = wp * w_1_11
838 750578400 : w10_00 = wp1 * w_0_00
839 750578400 : w10_01 = wp1 * w_0_01
840 750578400 : w10_10 = wp1 * w_0_10
841 750578400 : w10_11 = wp1 * w_0_11
842 750578400 : w11_00 = wp1 * w_1_00
843 750578400 : w11_01 = wp1 * w_1_01
844 750578400 : w11_10 = wp1 * w_1_10
845 750578400 : w11_11 = wp1 * w_1_11
846 : !
847 : ! Asymptotic value of absorptivity as U->infinity
848 : !
849 : fa = fat(1,ib) + &
850 : fat(2,ib) * te1 + &
851 : fat(3,ib) * te2 + &
852 : fat(4,ib) * te3 + &
853 : fat(5,ib) * te4 + &
854 750578400 : fat(6,ib) * te5
855 :
856 : a_star = &
857 3752892000 : ah2onw(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
858 750578400 : ah2onw(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
859 750578400 : ah2onw(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
860 : ah2onw(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
861 750578400 : ah2onw(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
862 : ah2onw(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
863 : ah2onw(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
864 : ah2onw(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
865 750578400 : ah2onw(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
866 : ah2onw(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
867 : ah2onw(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
868 : ah2onw(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
869 : ah2onw(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
870 : ah2onw(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
871 : ah2onw(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
872 : ah2onw(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
873 750578400 : ah2onw(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
874 : ah2onw(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
875 : ah2onw(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
876 : ah2onw(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
877 : ah2onw(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
878 : ah2onw(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
879 : ah2onw(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
880 : ah2onw(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
881 : ah2onw(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
882 : ah2onw(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
883 : ah2onw(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
884 : ah2onw(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
885 : ah2onw(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
886 : ah2onw(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
887 : ah2onw(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
888 7505784000 : ah2onw(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
889 : abso(i,ib) = min(max(fa * (1.0_r8 - (1.0_r8 - a_star) * &
890 : aer_trn_ttl(i,k1,k2,ib)), &
891 750578400 : 0.0_r8), 1.0_r8)
892 : !
893 : ! Invoke linear limit for scaling wrt u below min_u_h2o
894 : !
895 750578400 : if (uvar < min_u_h2o) then
896 7687876 : uscl = uvar / min_u_h2o
897 7687876 : abso(i,ib) = abso(i,ib) * uscl
898 : endif
899 :
900 : !
901 : ! Band-dependent indices for window
902 : !
903 750578400 : ib = 2
904 :
905 750578400 : uvar = ub(ib) * fdif
906 750578400 : log_u = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
907 750578400 : dvar = (log_u - min_lu_h2o) / dlu_h2o
908 750578400 : iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
909 750578400 : iu1 = iu + 1
910 750578400 : wu = dvar - floor(dvar)
911 750578400 : wu1 = 1.0_r8 - wu
912 :
913 750578400 : log_p = min(log10(max(pnewb(ib), min_p_h2o)), max_lp_h2o)
914 750578400 : dvar = (log_p - min_lp_h2o) / dlp_h2o
915 750578400 : ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
916 750578400 : ip1 = ip + 1
917 750578400 : wp = dvar - floor(dvar)
918 750578400 : wp1 = 1.0_r8 - wp
919 :
920 750578400 : w00_00 = wp * w_0_00
921 750578400 : w00_01 = wp * w_0_01
922 750578400 : w00_10 = wp * w_0_10
923 750578400 : w00_11 = wp * w_0_11
924 750578400 : w01_00 = wp * w_1_00
925 750578400 : w01_01 = wp * w_1_01
926 750578400 : w01_10 = wp * w_1_10
927 750578400 : w01_11 = wp * w_1_11
928 750578400 : w10_00 = wp1 * w_0_00
929 750578400 : w10_01 = wp1 * w_0_01
930 750578400 : w10_10 = wp1 * w_0_10
931 750578400 : w10_11 = wp1 * w_0_11
932 750578400 : w11_00 = wp1 * w_1_00
933 750578400 : w11_01 = wp1 * w_1_01
934 750578400 : w11_10 = wp1 * w_1_10
935 750578400 : w11_11 = wp1 * w_1_11
936 :
937 750578400 : log_uc = min(log10(max(uch2o * fdif, min_u_h2o)), max_lu_h2o)
938 750578400 : dvar = (log_uc - min_lu_h2o) / dlu_h2o
939 750578400 : iuc = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
940 750578400 : iuc1 = iuc + 1
941 750578400 : wuc = dvar - floor(dvar)
942 750578400 : wuc1 = 1.0_r8 - wuc
943 : !
944 : ! Asymptotic value of absorptivity as U->infinity
945 : !
946 : fa = fat(1,ib) + &
947 : fat(2,ib) * te1 + &
948 : fat(3,ib) * te2 + &
949 : fat(4,ib) * te3 + &
950 : fat(5,ib) * te4 + &
951 750578400 : fat(6,ib) * te5
952 :
953 : l_star = &
954 1501156800 : ln_ah2ow(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
955 : ln_ah2ow(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
956 : ln_ah2ow(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
957 : ln_ah2ow(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
958 750578400 : ln_ah2ow(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
959 : ln_ah2ow(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
960 : ln_ah2ow(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
961 : ln_ah2ow(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
962 : ln_ah2ow(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
963 : ln_ah2ow(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
964 : ln_ah2ow(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
965 : ln_ah2ow(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
966 : ln_ah2ow(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
967 : ln_ah2ow(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
968 : ln_ah2ow(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
969 : ln_ah2ow(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
970 750578400 : ln_ah2ow(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
971 : ln_ah2ow(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
972 : ln_ah2ow(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
973 : ln_ah2ow(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
974 : ln_ah2ow(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
975 : ln_ah2ow(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
976 : ln_ah2ow(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
977 : ln_ah2ow(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
978 : ln_ah2ow(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
979 : ln_ah2ow(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
980 : ln_ah2ow(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
981 : ln_ah2ow(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
982 : ln_ah2ow(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
983 : ln_ah2ow(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
984 : ln_ah2ow(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
985 3002313600 : ln_ah2ow(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
986 :
987 : c_star = &
988 750578400 : cn_ah2ow(ip , itp , iuc , ite , irh ) * w11_11 * wuc1 + &
989 : cn_ah2ow(ip , itp , iuc , ite , irh1) * w11_10 * wuc1 + &
990 : cn_ah2ow(ip , itp , iuc , ite1, irh ) * w11_01 * wuc1 + &
991 : cn_ah2ow(ip , itp , iuc , ite1, irh1) * w11_00 * wuc1 + &
992 750578400 : cn_ah2ow(ip , itp , iuc1, ite , irh ) * w11_11 * wuc + &
993 : cn_ah2ow(ip , itp , iuc1, ite , irh1) * w11_10 * wuc + &
994 : cn_ah2ow(ip , itp , iuc1, ite1, irh ) * w11_01 * wuc + &
995 : cn_ah2ow(ip , itp , iuc1, ite1, irh1) * w11_00 * wuc + &
996 : cn_ah2ow(ip , itp1, iuc , ite , irh ) * w10_11 * wuc1 + &
997 : cn_ah2ow(ip , itp1, iuc , ite , irh1) * w10_10 * wuc1 + &
998 : cn_ah2ow(ip , itp1, iuc , ite1, irh ) * w10_01 * wuc1 + &
999 : cn_ah2ow(ip , itp1, iuc , ite1, irh1) * w10_00 * wuc1 + &
1000 : cn_ah2ow(ip , itp1, iuc1, ite , irh ) * w10_11 * wuc + &
1001 : cn_ah2ow(ip , itp1, iuc1, ite , irh1) * w10_10 * wuc + &
1002 : cn_ah2ow(ip , itp1, iuc1, ite1, irh ) * w10_01 * wuc + &
1003 : cn_ah2ow(ip , itp1, iuc1, ite1, irh1) * w10_00 * wuc + &
1004 : cn_ah2ow(ip1, itp , iuc , ite , irh ) * w01_11 * wuc1 + &
1005 : cn_ah2ow(ip1, itp , iuc , ite , irh1) * w01_10 * wuc1 + &
1006 : cn_ah2ow(ip1, itp , iuc , ite1, irh ) * w01_01 * wuc1 + &
1007 : cn_ah2ow(ip1, itp , iuc , ite1, irh1) * w01_00 * wuc1 + &
1008 : cn_ah2ow(ip1, itp , iuc1, ite , irh ) * w01_11 * wuc + &
1009 : cn_ah2ow(ip1, itp , iuc1, ite , irh1) * w01_10 * wuc + &
1010 : cn_ah2ow(ip1, itp , iuc1, ite1, irh ) * w01_01 * wuc + &
1011 : cn_ah2ow(ip1, itp , iuc1, ite1, irh1) * w01_00 * wuc + &
1012 : cn_ah2ow(ip1, itp1, iuc , ite , irh ) * w00_11 * wuc1 + &
1013 : cn_ah2ow(ip1, itp1, iuc , ite , irh1) * w00_10 * wuc1 + &
1014 : cn_ah2ow(ip1, itp1, iuc , ite1, irh ) * w00_01 * wuc1 + &
1015 : cn_ah2ow(ip1, itp1, iuc , ite1, irh1) * w00_00 * wuc1 + &
1016 : cn_ah2ow(ip1, itp1, iuc1, ite , irh ) * w00_11 * wuc + &
1017 : cn_ah2ow(ip1, itp1, iuc1, ite , irh1) * w00_10 * wuc + &
1018 : cn_ah2ow(ip1, itp1, iuc1, ite1, irh ) * w00_01 * wuc + &
1019 1501156800 : cn_ah2ow(ip1, itp1, iuc1, ite1, irh1) * w00_00 * wuc
1020 : abso(i,ib) = min(max(fa * (1.0_r8 - l_star * c_star * &
1021 : aer_trn_ttl(i,k1,k2,ib)), &
1022 750578400 : 0.0_r8), 1.0_r8)
1023 : !
1024 : ! Invoke linear limit for scaling wrt u below min_u_h2o
1025 : !
1026 1548971424 : if (uvar < min_u_h2o) then
1027 7702230 : uscl = uvar / min_u_h2o
1028 7702230 : abso(i,ib) = abso(i,ib) * uscl
1029 : endif
1030 :
1031 : end do
1032 : !
1033 : ! Line transmission in 800-1000 and 1000-1200 cm-1 intervals
1034 : !
1035 798393024 : do i=1,ncol
1036 750578400 : term7(i,1) = coefj(1,1) + coefj(2,1)*dty(i)*(1._r8 + c16*dty(i))
1037 750578400 : term8(i,1) = coefk(1,1) + coefk(2,1)*dty(i)*(1._r8 + c17*dty(i))
1038 750578400 : term7(i,2) = coefj(1,2) + coefj(2,2)*dty(i)*(1._r8 + c26*dty(i))
1039 798393024 : term8(i,2) = coefk(1,2) + coefk(2,2)*dty(i)*(1._r8 + c27*dty(i))
1040 : end do
1041 : !
1042 : ! 500 - 800 cm-1 h2o rotation band overlap with co2
1043 : !
1044 798393024 : do i=1,ncol
1045 750578400 : k21 = term7(i,1) + term8(i,1)/ &
1046 750578400 : (1._r8 + (c30 + c31*(dty(i)-10._r8)*(dty(i)-10._r8))*sqrtu(i))
1047 : k22 = term7(i,2) + term8(i,2)/ &
1048 750578400 : (1._r8 + (c28 + c29*(dty(i)-10._r8))*sqrtu(i))
1049 750578400 : tr1 = exp(-(k21*(sqrtu(i) + fc1*fwku(i))))
1050 750578400 : tr2 = exp(-(k22*(sqrtu(i) + fc1*fwku(i))))
1051 750578400 : tr1=tr1*aer_trn_ttl(i,k1,k2,idx_LW_0650_0800)
1052 : ! ! H2O line+STRAER trn 650--800 cm-1
1053 750578400 : tr2=tr2*aer_trn_ttl(i,k1,k2,idx_LW_0500_0650)
1054 : ! ! H2O line+STRAER trn 500--650 cm-1
1055 750578400 : tr5 = exp(-((coefh(1,3) + coefh(2,3)*dtx(i))*uc1(i)))
1056 750578400 : tr6 = exp(-((coefh(1,4) + coefh(2,4)*dtx(i))*uc1(i)))
1057 750578400 : tr9(i) = tr1*tr5
1058 750578400 : tr10(i) = tr2*tr6
1059 750578400 : th2o(i) = tr10(i)
1060 798393024 : trab2(i) = 0.65_r8*tr9(i) + 0.35_r8*tr10(i)
1061 : end do
1062 47814624 : if (k2 < k1) then
1063 399196512 : do i=1,ncol
1064 399196512 : to3h2o(i) = h2otr(i,k1)/h2otr(i,k2)
1065 : end do
1066 : else
1067 399196512 : do i=1,ncol
1068 399196512 : to3h2o(i) = h2otr(i,k2)/h2otr(i,k1)
1069 : end do
1070 : end if
1071 : !
1072 : ! abso(i,3) o3 9.6 micrometer band (nu3 and nu1 bands)
1073 : !
1074 798393024 : do i=1,ncol
1075 750578400 : dpnm(i) = pnm(i,k1) - pnm(i,k2)
1076 750578400 : to3co2(i) = (pnm(i,k1)*co2t(i,k1) - pnm(i,k2)*co2t(i,k2))/dpnm(i)
1077 750578400 : te = (to3co2(i)*r293)**.7_r8
1078 750578400 : dplos = plos(i,k1) - plos(i,k2)
1079 798393024 : if (dplos == 0._r8) then
1080 0 : abso(i,3) = 0._r8
1081 0 : to3(i) = 1._r8
1082 0 : write(iulog,*) 'radiation ozone error ',k1,k2,plos(i,k1)
1083 : else
1084 750578400 : dplol = plol(i,k1) - plol(i,k2)
1085 750578400 : u1 = 18.29_r8*abs(dplos)/te
1086 750578400 : u2 = .5649_r8*abs(dplos)/te
1087 750578400 : rphat = dplol/dplos
1088 750578400 : tlocal = tint(i,k2)
1089 750578400 : tcrfac = sqrt(tlocal*r250)*te
1090 750578400 : beta = r3205*(rphat + dpfo3*tcrfac)
1091 750578400 : realnu = te/beta
1092 750578400 : tmp1 = u1/sqrt(4._r8 + u1*(1._r8 + realnu))
1093 750578400 : tmp2 = u2/sqrt(4._r8 + u2*(1._r8 + realnu))
1094 750578400 : o3bndi = 74._r8*te*log(1._r8 + tmp1 + tmp2)
1095 750578400 : abso(i,3) = o3bndi*to3h2o(i)*dbvtit(i,k2)
1096 750578400 : to3(i) = 1.0_r8/(1._r8 + 0.1_r8*tmp1 + 0.1_r8*tmp2)
1097 : endif
1098 : end do
1099 : !
1100 : ! abso(i,4) co2 15 micrometer band system
1101 : !
1102 798393024 : do i=1,ncol
1103 750578400 : sqwp = sqrt(abs(plco2(i,k1) - plco2(i,k2)))
1104 750578400 : et = exp(-480._r8/to3co2(i))
1105 750578400 : sqti(i) = sqrt(to3co2(i))
1106 750578400 : rsqti = 1._r8/sqti(i)
1107 750578400 : et2 = et*et
1108 750578400 : et4 = et2*et2
1109 750578400 : omet = 1._r8 - 1.5_r8*et2
1110 750578400 : f1co2 = 899.70_r8*omet*(1._r8 + 1.94774_r8*et + 4.73486_r8*et2)*rsqti
1111 750578400 : f1sqwp(i) = f1co2*sqwp
1112 750578400 : t1co2(i) = 1._r8/(1._r8 + (245.18_r8*omet*sqwp*rsqti))
1113 750578400 : oneme = 1._r8 - et2
1114 750578400 : alphat = oneme**3*rsqti
1115 750578400 : pi = abs(dpnm(i))
1116 750578400 : wco2 = 2.5221_r8*co2vmr(i)*pi*rga
1117 750578400 : u7(i) = 4.9411e4_r8*alphat*et2*wco2
1118 750578400 : u8 = 3.9744e4_r8*alphat*et4*wco2
1119 750578400 : u9 = 1.0447e5_r8*alphat*et4*et2*wco2
1120 750578400 : u13 = 2.8388e3_r8*alphat*et4*wco2
1121 750578400 : tpath = to3co2(i)
1122 750578400 : tlocal = tint(i,k2)
1123 750578400 : tcrfac = sqrt(tlocal*r250*tpath*r300)
1124 750578400 : posqt = ((pnm(i,k2) + pnm(i,k1))*r2sslp + dpfco2*tcrfac)*rsqti
1125 750578400 : rbeta7(i) = 1._r8/(5.3228_r8*posqt)
1126 750578400 : rbeta8 = 1._r8/(10.6576_r8*posqt)
1127 750578400 : rbeta9 = rbeta7(i)
1128 750578400 : rbeta13 = rbeta9
1129 : f2co2(i) = (u7(i)/sqrt(4._r8 + u7(i)*(1._r8 + rbeta7(i)))) + &
1130 : (u8 /sqrt(4._r8 + u8*(1._r8 + rbeta8))) + &
1131 750578400 : (u9 /sqrt(4._r8 + u9*(1._r8 + rbeta9)))
1132 798393024 : f3co2(i) = u13/sqrt(4._r8 + u13*(1._r8 + rbeta13))
1133 : end do
1134 47814624 : if (k2 >= k1) then
1135 399196512 : do i=1,ncol
1136 399196512 : sqti(i) = sqrt(tlayr(i,k2))
1137 : end do
1138 : end if
1139 : !
1140 798393024 : do i=1,ncol
1141 750578400 : tmp1 = log(1._r8 + f1sqwp(i))
1142 750578400 : tmp2 = log(1._r8 + f2co2(i))
1143 750578400 : tmp3 = log(1._r8 + f3co2(i))
1144 750578400 : absbnd = (tmp1 + 2._r8*t1co2(i)*tmp2 + 2._r8*tmp3)*sqti(i)
1145 750578400 : abso(i,4) = trab2(i)*co2em(i,k2)*absbnd
1146 798393024 : tco2(i) = 1._r8/(1.0_r8+10.0_r8*(u7(i)/sqrt(4._r8 + u7(i)*(1._r8 + rbeta7(i)))))
1147 : end do
1148 : !
1149 : ! Calculate absorptivity due to trace gases, abstrc
1150 : !
1151 : call trcab(ncol , &
1152 : k1 ,k2 ,ucfc11 ,ucfc12 ,un2o0 , &
1153 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
1154 : uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
1155 : bch4 ,to3co2 ,pnm ,dw ,pnew , &
1156 : s2c ,uptype ,u ,abplnk1 ,tco2 , &
1157 : th2o ,to3 ,abstrc , &
1158 47814624 : aer_trn_ttl)
1159 : !
1160 : ! Sum total absorptivity
1161 : !
1162 800232048 : do i=1,ncol
1163 2251735200 : abstot(i,k1,k2) = abso(i,1) + abso(i,2) + &
1164 3051967248 : abso(i,3) + abso(i,4) + abstrc(i)
1165 : end do
1166 : end do ! do k2 =
1167 : end do ! do k1 =
1168 : !
1169 : ! Adjacent layer absorptivity:
1170 : !
1171 : ! abso(i,1) 0 - 800 cm-1 h2o rotation band
1172 : ! abso(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
1173 : ! abso(i,2) 800 - 1200 cm-1 h2o window
1174 : !
1175 : ! Separation between rotation and vibration-rotation dropped, so
1176 : ! only 2 slots needed for H2O absorptivity
1177 : !
1178 : ! 500-800 cm^-1 H2o continuum/line overlap already included
1179 : ! in abso(i,1). This used to be in abso(i,4)
1180 : !
1181 : ! abso(i,3) o3 9.6 micrometer band (nu3 and nu1 bands)
1182 : ! abso(i,4) co2 15 micrometer band system
1183 : !
1184 : ! Nearest layer level loop
1185 : !
1186 1839024 : do k2=pver,ntoplw,-1
1187 29570112 : do i=1,ncol
1188 27799200 : tbar(i,1) = 0.5_r8*(tint(i,k2+1) + tlayr(i,k2+1))
1189 27799200 : emm(i,1) = 0.5_r8*(co2em(i,k2+1) + co2eml(i,k2))
1190 27799200 : tbar(i,2) = 0.5_r8*(tlayr(i,k2+1) + tint(i,k2))
1191 27799200 : emm(i,2) = 0.5_r8*(co2em(i,k2) + co2eml(i,k2))
1192 27799200 : tbar(i,3) = 0.5_r8*(tbar(i,2) + tbar(i,1))
1193 27799200 : emm(i,3) = emm(i,1)
1194 27799200 : tbar(i,4) = tbar(i,3)
1195 27799200 : emm(i,4) = emm(i,2)
1196 27799200 : o3emm(i,1) = 0.5_r8*(dbvtit(i,k2+1) + dbvtly(i,k2))
1197 27799200 : o3emm(i,2) = 0.5_r8*(dbvtit(i,k2) + dbvtly(i,k2))
1198 27799200 : o3emm(i,3) = o3emm(i,1)
1199 27799200 : o3emm(i,4) = o3emm(i,2)
1200 27799200 : temh2o(i,1) = tbar(i,1)
1201 27799200 : temh2o(i,2) = tbar(i,2)
1202 27799200 : temh2o(i,3) = tbar(i,1)
1203 27799200 : temh2o(i,4) = tbar(i,2)
1204 29570112 : dpnm(i) = pnm(i,k2+1) - pnm(i,k2)
1205 : end do
1206 : !
1207 : ! Weighted Planck functions for trace gases
1208 : !
1209 26563680 : do wvl = 1,14
1210 415752480 : do i = 1,ncol
1211 389188800 : bplnk(wvl,i,1) = 0.5_r8*(abplnk1(wvl,i,k2+1) + abplnk2(wvl,i,k2))
1212 389188800 : bplnk(wvl,i,2) = 0.5_r8*(abplnk1(wvl,i,k2) + abplnk2(wvl,i,k2))
1213 389188800 : bplnk(wvl,i,3) = bplnk(wvl,i,1)
1214 413981568 : bplnk(wvl,i,4) = bplnk(wvl,i,2)
1215 : end do
1216 : end do
1217 :
1218 29570112 : do i=1,ncol
1219 27799200 : rdpnmsq = 1._r8/(pnmsq(i,k2+1) - pnmsq(i,k2))
1220 27799200 : rdpnm = 1._r8/dpnm(i)
1221 27799200 : p1 = .5_r8*(pbr(i,k2) + pnm(i,k2+1))
1222 27799200 : p2 = .5_r8*(pbr(i,k2) + pnm(i,k2 ))
1223 27799200 : uinpl(i,1) = (pnmsq(i,k2+1) - p1**2)*rdpnmsq
1224 27799200 : uinpl(i,2) = -(pnmsq(i,k2 ) - p2**2)*rdpnmsq
1225 27799200 : uinpl(i,3) = -(pnmsq(i,k2 ) - p1**2)*rdpnmsq
1226 27799200 : uinpl(i,4) = (pnmsq(i,k2+1) - p2**2)*rdpnmsq
1227 27799200 : winpl(i,1) = (.5_r8*( pnm(i,k2+1) - pbr(i,k2)))*rdpnm
1228 27799200 : winpl(i,2) = (.5_r8*(-pnm(i,k2 ) + pbr(i,k2)))*rdpnm
1229 27799200 : winpl(i,3) = (.5_r8*( pnm(i,k2+1) + pbr(i,k2)) - pnm(i,k2 ))*rdpnm
1230 27799200 : winpl(i,4) = (.5_r8*(-pnm(i,k2 ) - pbr(i,k2)) + pnm(i,k2+1))*rdpnm
1231 27799200 : tmp1 = 1._r8/(piln(i,k2+1) - piln(i,k2))
1232 27799200 : tmp2 = piln(i,k2+1) - pmln(i,k2)
1233 27799200 : tmp3 = piln(i,k2 ) - pmln(i,k2)
1234 27799200 : zinpl(i,1) = (.5_r8*tmp2 )*tmp1
1235 27799200 : zinpl(i,2) = ( - .5_r8*tmp3)*tmp1
1236 27799200 : zinpl(i,3) = (.5_r8*tmp2 - tmp3)*tmp1
1237 27799200 : zinpl(i,4) = ( tmp2 - .5_r8*tmp3)*tmp1
1238 27799200 : pinpl(i,1) = 0.5_r8*(p1 + pnm(i,k2+1))
1239 27799200 : pinpl(i,2) = 0.5_r8*(p2 + pnm(i,k2 ))
1240 27799200 : pinpl(i,3) = 0.5_r8*(p1 + pnm(i,k2 ))
1241 29570112 : pinpl(i,4) = 0.5_r8*(p2 + pnm(i,k2+1))
1242 : end do
1243 8922672 : do kn=1,4
1244 118280448 : do i=1,ncol
1245 111196800 : u(i) = uinpl(i,kn)*abs(plh2o(i,k2) - plh2o(i,k2+1))
1246 111196800 : sqrtu(i) = sqrt(u(i))
1247 111196800 : dw(i) = abs(w(i,k2) - w(i,k2+1))
1248 111196800 : pnew(i) = u(i)/(winpl(i,kn)*dw(i))
1249 111196800 : pnew_mks = pnew(i) * sslp_mks
1250 111196800 : t_p = min(max(tbar(i,kn), min_tp_h2o), max_tp_h2o)
1251 :
1252 111196800 : call qsat_water(t_p, pnew_mks, esx, qsx)
1253 :
1254 111196800 : q_path = dw(i) / ABS(dpnm(i)) / rga
1255 :
1256 111196800 : ds2c = abs(s2c(i,k2) - s2c(i,k2+1))
1257 111196800 : uc1(i) = uinpl(i,kn)*ds2c
1258 111196800 : pch2o = uc1(i)
1259 111196800 : uc1(i) = (uc1(i) + 1.7e-3_r8*u(i))*(1._r8 + 2._r8*uc1(i))/(1._r8 + 15._r8*uc1(i))
1260 111196800 : dtx(i) = temh2o(i,kn) - 250._r8
1261 111196800 : dty(i) = tbar(i,kn) - 250._r8
1262 :
1263 111196800 : fwk(i) = fwcoef + fwc1/(1._r8 + fwc2*u(i))
1264 111196800 : fwku(i) = fwk(i)*u(i)
1265 :
1266 : aer_trn_ngh(i, 1:nlwbands)= &
1267 889574400 : exp(-fdif * uinpl(i,kn) * odap_aer(i, k2, 1:nlwbands ) )
1268 :
1269 : !
1270 : ! Define variables for C/H/E (now C/LT/E) fit
1271 : !
1272 : ! abso(i,1) 0 - 800 cm-1 h2o rotation band
1273 : ! abso(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
1274 : ! abso(i,2) 800 - 1200 cm-1 h2o window
1275 : !
1276 : ! Separation between rotation and vibration-rotation dropped, so
1277 : ! only 2 slots needed for H2O absorptivity
1278 : !
1279 : ! Notation:
1280 : ! U = integral (P/P_0 dW)
1281 : ! P = atmospheric pressure
1282 : ! P_0 = reference atmospheric pressure
1283 : ! W = precipitable water path
1284 : ! T_e = emission temperature
1285 : ! T_p = path temperature
1286 : ! RH = path relative humidity
1287 : !
1288 : !
1289 : ! Terms for asymptotic value of emissivity
1290 : !
1291 111196800 : te1 = temh2o(i,kn)
1292 111196800 : te2 = te1 * te1
1293 111196800 : te3 = te2 * te1
1294 111196800 : te4 = te3 * te1
1295 111196800 : te5 = te4 * te1
1296 :
1297 : !
1298 : ! Indices for lines and continuum tables
1299 : ! Note: because we are dealing with the nearest layer,
1300 : ! the Hulst-Curtis-Godson corrections
1301 : ! for inhomogeneous paths are not applied.
1302 : !
1303 111196800 : uvar = u(i)*fdif
1304 111196800 : log_u = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
1305 111196800 : dvar = (log_u - min_lu_h2o) / dlu_h2o
1306 111196800 : iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
1307 111196800 : iu1 = iu + 1
1308 111196800 : wu = dvar - floor(dvar)
1309 111196800 : wu1 = 1.0_r8 - wu
1310 :
1311 111196800 : log_p = min(log10(max(pnew(i), min_p_h2o)), max_lp_h2o)
1312 111196800 : dvar = (log_p - min_lp_h2o) / dlp_h2o
1313 111196800 : ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
1314 111196800 : ip1 = ip + 1
1315 111196800 : wp = dvar - floor(dvar)
1316 111196800 : wp1 = 1.0_r8 - wp
1317 :
1318 111196800 : dvar = (t_p - min_tp_h2o) / dtp_h2o
1319 111196800 : itp = min(max(int(aint(dvar,r8)) + 1, 1), n_tp - 1)
1320 111196800 : itp1 = itp + 1
1321 111196800 : wtp = dvar - floor(dvar)
1322 111196800 : wtp1 = 1.0_r8 - wtp
1323 :
1324 111196800 : t_e = min(max(temh2o(i,kn)-t_p,min_te_h2o),max_te_h2o)
1325 111196800 : dvar = (t_e - min_te_h2o) / dte_h2o
1326 111196800 : ite = min(max(int(aint(dvar,r8)) + 1, 1), n_te - 1)
1327 111196800 : ite1 = ite + 1
1328 111196800 : wte = dvar - floor(dvar)
1329 111196800 : wte1 = 1.0_r8 - wte
1330 :
1331 111196800 : rh_path = min(max(q_path / qsx, min_rh_h2o), max_rh_h2o)
1332 111196800 : dvar = (rh_path - min_rh_h2o) / drh_h2o
1333 111196800 : irh = min(max(int(aint(dvar,r8)) + 1, 1), n_rh - 1)
1334 111196800 : irh1 = irh + 1
1335 111196800 : wrh = dvar - floor(dvar)
1336 111196800 : wrh1 = 1.0_r8 - wrh
1337 :
1338 111196800 : w_0_0_ = wtp * wte
1339 111196800 : w_0_1_ = wtp * wte1
1340 111196800 : w_1_0_ = wtp1 * wte
1341 111196800 : w_1_1_ = wtp1 * wte1
1342 :
1343 111196800 : w_0_00 = w_0_0_ * wrh
1344 111196800 : w_0_01 = w_0_0_ * wrh1
1345 111196800 : w_0_10 = w_0_1_ * wrh
1346 111196800 : w_0_11 = w_0_1_ * wrh1
1347 111196800 : w_1_00 = w_1_0_ * wrh
1348 111196800 : w_1_01 = w_1_0_ * wrh1
1349 111196800 : w_1_10 = w_1_1_ * wrh
1350 111196800 : w_1_11 = w_1_1_ * wrh1
1351 :
1352 111196800 : w00_00 = wp * w_0_00
1353 111196800 : w00_01 = wp * w_0_01
1354 111196800 : w00_10 = wp * w_0_10
1355 111196800 : w00_11 = wp * w_0_11
1356 111196800 : w01_00 = wp * w_1_00
1357 111196800 : w01_01 = wp * w_1_01
1358 111196800 : w01_10 = wp * w_1_10
1359 111196800 : w01_11 = wp * w_1_11
1360 111196800 : w10_00 = wp1 * w_0_00
1361 111196800 : w10_01 = wp1 * w_0_01
1362 111196800 : w10_10 = wp1 * w_0_10
1363 111196800 : w10_11 = wp1 * w_0_11
1364 111196800 : w11_00 = wp1 * w_1_00
1365 111196800 : w11_01 = wp1 * w_1_01
1366 111196800 : w11_10 = wp1 * w_1_10
1367 111196800 : w11_11 = wp1 * w_1_11
1368 :
1369 : !
1370 : ! Non-window absorptivity
1371 : !
1372 111196800 : ib = 1
1373 :
1374 : fa = fat(1,ib) + &
1375 : fat(2,ib) * te1 + &
1376 : fat(3,ib) * te2 + &
1377 : fat(4,ib) * te3 + &
1378 : fat(5,ib) * te4 + &
1379 111196800 : fat(6,ib) * te5
1380 :
1381 : a_star = &
1382 555984000 : ah2onw(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
1383 111196800 : ah2onw(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
1384 111196800 : ah2onw(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
1385 : ah2onw(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
1386 111196800 : ah2onw(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
1387 : ah2onw(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
1388 : ah2onw(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
1389 : ah2onw(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
1390 111196800 : ah2onw(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
1391 : ah2onw(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
1392 : ah2onw(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
1393 : ah2onw(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
1394 : ah2onw(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
1395 : ah2onw(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
1396 : ah2onw(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
1397 : ah2onw(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
1398 111196800 : ah2onw(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
1399 : ah2onw(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
1400 : ah2onw(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
1401 : ah2onw(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
1402 : ah2onw(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
1403 : ah2onw(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
1404 : ah2onw(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
1405 : ah2onw(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
1406 : ah2onw(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
1407 : ah2onw(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
1408 : ah2onw(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
1409 : ah2onw(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
1410 : ah2onw(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
1411 : ah2onw(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
1412 : ah2onw(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
1413 1111968000 : ah2onw(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
1414 :
1415 : abso(i,ib) = min(max(fa * (1.0_r8 - (1.0_r8 - a_star) * &
1416 : aer_trn_ngh(i,ib)), &
1417 111196800 : 0.0_r8), 1.0_r8)
1418 :
1419 : !
1420 : ! Invoke linear limit for scaling wrt u below min_u_h2o
1421 : !
1422 111196800 : if (uvar < min_u_h2o) then
1423 10520652 : uscl = uvar / min_u_h2o
1424 10520652 : abso(i,ib) = abso(i,ib) * uscl
1425 : endif
1426 :
1427 : !
1428 : ! Window absorptivity
1429 : !
1430 111196800 : ib = 2
1431 :
1432 : fa = fat(1,ib) + &
1433 : fat(2,ib) * te1 + &
1434 : fat(3,ib) * te2 + &
1435 : fat(4,ib) * te3 + &
1436 : fat(5,ib) * te4 + &
1437 111196800 : fat(6,ib) * te5
1438 :
1439 : a_star = &
1440 : ah2ow(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
1441 : ah2ow(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
1442 : ah2ow(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
1443 : ah2ow(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
1444 : ah2ow(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
1445 : ah2ow(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
1446 : ah2ow(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
1447 : ah2ow(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
1448 : ah2ow(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
1449 : ah2ow(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
1450 : ah2ow(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
1451 : ah2ow(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
1452 : ah2ow(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
1453 : ah2ow(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
1454 : ah2ow(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
1455 : ah2ow(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
1456 : ah2ow(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
1457 : ah2ow(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
1458 : ah2ow(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
1459 : ah2ow(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
1460 : ah2ow(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
1461 : ah2ow(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
1462 : ah2ow(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
1463 : ah2ow(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
1464 : ah2ow(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
1465 : ah2ow(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
1466 : ah2ow(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
1467 : ah2ow(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
1468 : ah2ow(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
1469 : ah2ow(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
1470 : ah2ow(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
1471 111196800 : ah2ow(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
1472 :
1473 : abso(i,ib) = min(max(fa * (1.0_r8 - (1.0_r8 - a_star) * &
1474 : aer_trn_ngh(i,ib)), &
1475 111196800 : 0.0_r8), 1.0_r8)
1476 :
1477 : !
1478 : ! Invoke linear limit for scaling wrt u below min_u_h2o
1479 : !
1480 229477248 : if (uvar < min_u_h2o) then
1481 10520652 : uscl = uvar / min_u_h2o
1482 10520652 : abso(i,ib) = abso(i,ib) * uscl
1483 : endif
1484 :
1485 : end do
1486 : !
1487 : ! Line transmission in 800-1000 and 1000-1200 cm-1 intervals
1488 : !
1489 118280448 : do i=1,ncol
1490 111196800 : term7(i,1) = coefj(1,1) + coefj(2,1)*dty(i)*(1._r8 + c16*dty(i))
1491 111196800 : term8(i,1) = coefk(1,1) + coefk(2,1)*dty(i)*(1._r8 + c17*dty(i))
1492 111196800 : term7(i,2) = coefj(1,2) + coefj(2,2)*dty(i)*(1._r8 + c26*dty(i))
1493 118280448 : term8(i,2) = coefk(1,2) + coefk(2,2)*dty(i)*(1._r8 + c27*dty(i))
1494 : end do
1495 : !
1496 : ! 500 - 800 cm-1 h2o rotation band overlap with co2
1497 : !
1498 118280448 : do i=1,ncol
1499 111196800 : dtym10 = dty(i) - 10._r8
1500 111196800 : denom = 1._r8 + (c30 + c31*dtym10*dtym10)*sqrtu(i)
1501 111196800 : k21 = term7(i,1) + term8(i,1)/denom
1502 111196800 : denom = 1._r8 + (c28 + c29*dtym10 )*sqrtu(i)
1503 111196800 : k22 = term7(i,2) + term8(i,2)/denom
1504 111196800 : tr1 = exp(-(k21*(sqrtu(i) + fc1*fwku(i))))
1505 111196800 : tr2 = exp(-(k22*(sqrtu(i) + fc1*fwku(i))))
1506 111196800 : tr1=tr1*aer_trn_ngh(i,idx_LW_0650_0800)
1507 : ! ! H2O line+STRAER trn 650--800 cm-1
1508 111196800 : tr2=tr2*aer_trn_ngh(i,idx_LW_0500_0650)
1509 : ! ! H2O line+STRAER trn 500--650 cm-1
1510 111196800 : tr5 = exp(-((coefh(1,3) + coefh(2,3)*dtx(i))*uc1(i)))
1511 111196800 : tr6 = exp(-((coefh(1,4) + coefh(2,4)*dtx(i))*uc1(i)))
1512 111196800 : tr9(i) = tr1*tr5
1513 111196800 : tr10(i) = tr2*tr6
1514 111196800 : trab2(i)= 0.65_r8*tr9(i) + 0.35_r8*tr10(i)
1515 118280448 : th2o(i) = tr10(i)
1516 : end do
1517 : !
1518 : ! abso(i,3) o3 9.6 micrometer (nu3 and nu1 bands)
1519 : !
1520 118280448 : do i=1,ncol
1521 111196800 : te = (tbar(i,kn)*r293)**.7_r8
1522 111196800 : dplos = abs(plos(i,k2+1) - plos(i,k2))
1523 111196800 : u1 = zinpl(i,kn)*18.29_r8*dplos/te
1524 111196800 : u2 = zinpl(i,kn)*.5649_r8*dplos/te
1525 111196800 : tlocal = tbar(i,kn)
1526 111196800 : tcrfac = sqrt(tlocal*r250)*te
1527 111196800 : beta = r3205*(pinpl(i,kn)*rsslp + dpfo3*tcrfac)
1528 111196800 : realnu = te/beta
1529 111196800 : tmp1 = u1/sqrt(4._r8 + u1*(1._r8 + realnu))
1530 111196800 : tmp2 = u2/sqrt(4._r8 + u2*(1._r8 + realnu))
1531 111196800 : o3bndi = 74._r8*te*log(1._r8 + tmp1 + tmp2)
1532 111196800 : abso(i,3) = o3bndi*o3emm(i,kn)*(h2otr(i,k2+1)/h2otr(i,k2))
1533 118280448 : to3(i) = 1.0_r8/(1._r8 + 0.1_r8*tmp1 + 0.1_r8*tmp2)
1534 : end do
1535 : !
1536 : ! abso(i,4) co2 15 micrometer band system
1537 : !
1538 118280448 : do i=1,ncol
1539 111196800 : dplco2 = plco2(i,k2+1) - plco2(i,k2)
1540 111196800 : sqwp = sqrt(uinpl(i,kn)*dplco2)
1541 111196800 : et = exp(-480._r8/tbar(i,kn))
1542 111196800 : sqti(i) = sqrt(tbar(i,kn))
1543 111196800 : rsqti = 1._r8/sqti(i)
1544 111196800 : et2 = et*et
1545 111196800 : et4 = et2*et2
1546 111196800 : omet = (1._r8 - 1.5_r8*et2)
1547 111196800 : f1co2 = 899.70_r8*omet*(1._r8 + 1.94774_r8*et + 4.73486_r8*et2)*rsqti
1548 111196800 : f1sqwp(i)= f1co2*sqwp
1549 111196800 : t1co2(i) = 1._r8/(1._r8 + (245.18_r8*omet*sqwp*rsqti))
1550 111196800 : oneme = 1._r8 - et2
1551 111196800 : alphat = oneme**3*rsqti
1552 111196800 : pi = abs(dpnm(i))*winpl(i,kn)
1553 111196800 : wco2 = 2.5221_r8*co2vmr(i)*pi*rga
1554 111196800 : u7(i) = 4.9411e4_r8*alphat*et2*wco2
1555 111196800 : u8 = 3.9744e4_r8*alphat*et4*wco2
1556 111196800 : u9 = 1.0447e5_r8*alphat*et4*et2*wco2
1557 111196800 : u13 = 2.8388e3_r8*alphat*et4*wco2
1558 111196800 : tpath = tbar(i,kn)
1559 111196800 : tlocal = tbar(i,kn)
1560 111196800 : tcrfac = sqrt((tlocal*r250)*(tpath*r300))
1561 111196800 : posqt = (pinpl(i,kn)*rsslp + dpfco2*tcrfac)*rsqti
1562 111196800 : rbeta7(i)= 1._r8/(5.3228_r8*posqt)
1563 111196800 : rbeta8 = 1._r8/(10.6576_r8*posqt)
1564 111196800 : rbeta9 = rbeta7(i)
1565 111196800 : rbeta13 = rbeta9
1566 : f2co2(i) = u7(i)/sqrt(4._r8 + u7(i)*(1._r8 + rbeta7(i))) + &
1567 : u8 /sqrt(4._r8 + u8*(1._r8 + rbeta8)) + &
1568 111196800 : u9 /sqrt(4._r8 + u9*(1._r8 + rbeta9))
1569 111196800 : f3co2(i) = u13/sqrt(4._r8 + u13*(1._r8 + rbeta13))
1570 111196800 : tmp1 = log(1._r8 + f1sqwp(i))
1571 111196800 : tmp2 = log(1._r8 + f2co2(i))
1572 111196800 : tmp3 = log(1._r8 + f3co2(i))
1573 111196800 : absbnd = (tmp1 + 2._r8*t1co2(i)*tmp2 + 2._r8*tmp3)*sqti(i)
1574 111196800 : abso(i,4)= trab2(i)*emm(i,kn)*absbnd
1575 118280448 : tco2(i) = 1.0_r8/(1.0_r8+ 10.0_r8*u7(i)/sqrt(4._r8 + u7(i)*(1._r8 + rbeta7(i))))
1576 : end do ! do i =
1577 : !
1578 : ! Calculate trace gas absorptivity for nearest layer, abstrc
1579 : !
1580 : call trcabn(ncol , &
1581 : k2 ,kn ,ucfc11 ,ucfc12 ,un2o0 , &
1582 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
1583 : uco221 ,uco222 ,uco223 ,tbar ,bplnk , &
1584 : winpl ,pinpl ,tco2 ,th2o ,to3 , &
1585 : uptype ,dw ,s2c ,u ,pnew , &
1586 : abstrc ,uinpl , &
1587 7083648 : aer_trn_ngh)
1588 : !
1589 : ! Total next layer absorptivity:
1590 : !
1591 120051360 : do i=1,ncol
1592 333590400 : absnxt(i,k2,kn) = abso(i,1) + abso(i,2) + &
1593 451870848 : abso(i,3) + abso(i,4) + abstrc(i)
1594 : end do
1595 : end do ! do kn =
1596 : end do ! do k2 =
1597 :
1598 68112 : return
1599 : end subroutine radabs
1600 :
1601 : !====================================================================================
1602 :
1603 68112 : subroutine radems(lchnk ,ncol , &
1604 : s2c ,tcg ,w ,tplnke ,plh2o , &
1605 : pnm ,plco2 ,tint ,tint4 ,tlayr , &
1606 : tlayr4 ,plol ,plos ,ucfc11 ,ucfc12 , &
1607 : un2o0 ,un2o1 ,uch4 ,uco211 ,uco212 , &
1608 : uco213 ,uco221 ,uco222 ,uco223 ,uptype , &
1609 : bn2o0 ,bn2o1 ,bch4 ,co2em ,co2eml , &
1610 : co2t ,h2otr ,abplnk1 ,abplnk2 ,emstot , &
1611 : plh2ob ,wb , &
1612 : aer_trn_ttl, co2mmr)
1613 : !-----------------------------------------------------------------------
1614 : !
1615 : ! Purpose:
1616 : ! Compute emissivity for H2O, CO2, O3, CH4, N2O, CFC11 and CFC12
1617 : !
1618 : ! Method:
1619 : ! H2O .... Uses nonisothermal emissivity method for water vapor from
1620 : ! Ramanathan, V. and P.Downey, 1986: A Nonisothermal
1621 : ! Emissivity and Absorptivity Formulation for Water Vapor
1622 : ! Jouranl of Geophysical Research, vol. 91., D8, pp 8649-8666
1623 : !
1624 : ! Implementation updated by Collins,Hackney, and Edwards 2001
1625 : ! using line-by-line calculations based upon Hitran 1996 and
1626 : ! CKD 2.1 for absorptivity and emissivity
1627 : !
1628 : ! Implementation updated by Collins, Lee-Taylor, and Edwards (2003)
1629 : ! using line-by-line calculations based upon Hitran 2000 and
1630 : ! CKD 2.4 for absorptivity and emissivity
1631 : !
1632 : ! CO2 .... Uses absorptance parameterization of the 15 micro-meter
1633 : ! (500 - 800 cm-1) band system of Carbon Dioxide, from
1634 : ! Kiehl, J.T. and B.P.Briegleb, 1991: A New Parameterization
1635 : ! of the Absorptance Due to the 15 micro-meter Band System
1636 : ! of Carbon Dioxide Jouranl of Geophysical Research,
1637 : ! vol. 96., D5, pp 9013-9019. Also includes the effects
1638 : ! of the 9.4 and 10.4 micron bands of CO2.
1639 : !
1640 : ! O3 .... Uses absorptance parameterization of the 9.6 micro-meter
1641 : ! band system of ozone, from Ramanathan, V. and R. Dickinson,
1642 : ! 1979: The Role of stratospheric ozone in the zonal and
1643 : ! seasonal radiative energy balance of the earth-troposphere
1644 : ! system. Journal of the Atmospheric Sciences, Vol. 36,
1645 : ! pp 1084-1104
1646 : !
1647 : ! ch4 .... Uses a broad band model for the 7.7 micron band of methane.
1648 : !
1649 : ! n20 .... Uses a broad band model for the 7.8, 8.6 and 17.0 micron
1650 : ! bands of nitrous oxide
1651 : !
1652 : ! cfc11 ... Uses a quasi-linear model for the 9.2, 10.7, 11.8 and 12.5
1653 : ! micron bands of CFC11
1654 : !
1655 : ! cfc12 ... Uses a quasi-linear model for the 8.6, 9.1, 10.8 and 11.2
1656 : ! micron bands of CFC12
1657 : !
1658 : !
1659 : ! Computes individual emissivities, accounting for band overlap, and
1660 : ! sums to obtain the total.
1661 : !
1662 : ! Author: W. Collins (H2O emissivity) and J. Kiehl
1663 : !------------------------------Arguments--------------------------------
1664 : !
1665 : ! Input arguments
1666 : !
1667 : integer, intent(in) :: lchnk ! chunk identifier
1668 : integer, intent(in) :: ncol ! number of atmospheric columns
1669 :
1670 : real(r8), intent(in) :: s2c(pcols,pverp) ! H2o continuum path length
1671 : real(r8), intent(in) :: tcg(pcols,pverp) ! H2o-mass-wgted temp. (Curtis-Godson approx.)
1672 : real(r8), intent(in) :: w(pcols,pverp) ! H2o path length
1673 : real(r8), intent(in) :: tplnke(pcols) ! Layer planck temperature
1674 : real(r8), intent(in) :: plh2o(pcols,pverp) ! H2o prs wghted path length
1675 : real(r8), intent(in) :: pnm(pcols,pverp) ! Model interface pressure
1676 : real(r8), intent(in) :: plco2(pcols,pverp) ! Prs wghted path of co2
1677 : real(r8), intent(in) :: tint(pcols,pverp) ! Model interface temperatures
1678 : real(r8), intent(in) :: tint4(pcols,pverp) ! Tint to the 4th power
1679 : real(r8), intent(in) :: tlayr(pcols,pverp) ! K-1 model layer temperature
1680 : real(r8), intent(in) :: tlayr4(pcols,pverp) ! Tlayr to the 4th power
1681 : real(r8), intent(in) :: plol(pcols,pverp) ! Pressure wghtd ozone path
1682 : real(r8), intent(in) :: plos(pcols,pverp) ! Ozone path
1683 : real(r8), intent(in) :: plh2ob(nbands,pcols,pverp) ! Pressure weighted h2o path with
1684 : ! Hulst-Curtis-Godson temp. factor
1685 : ! for H2O bands
1686 : real(r8), intent(in) :: wb(nbands,pcols,pverp) ! H2o path length with
1687 : ! Hulst-Curtis-Godson temp. factor
1688 : ! for H2O bands
1689 :
1690 : real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands)
1691 : ! ! [fraction] Total strat. aerosol
1692 : ! ! transmission between interfaces k1 and k2
1693 :
1694 : !
1695 : ! Trace gas variables
1696 : !
1697 : real(r8), intent(in) :: co2mmr(pcols) ! co2 column mean mass mixing ratio
1698 : real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
1699 : real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
1700 : real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
1701 : real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
1702 : real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
1703 : real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
1704 : real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
1705 : real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
1706 : real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
1707 : real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
1708 : real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
1709 : real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
1710 : real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
1711 : real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
1712 : real(r8), intent(in) :: uptype(pcols,pverp) ! p-type continuum path length
1713 : !
1714 : ! Output arguments
1715 : !
1716 : real(r8), intent(out) :: emstot(pcols,pverp) ! Total emissivity
1717 : real(r8), intent(out) :: co2em(pcols,pverp) ! Layer co2 normalzd plnck funct drvtv
1718 : real(r8), intent(out) :: co2eml(pcols,pver) ! Intrfc co2 normalzd plnck func drvtv
1719 : real(r8), intent(out) :: co2t(pcols,pverp) ! Tmp and prs weighted path length
1720 : real(r8), intent(out) :: h2otr(pcols,pverp) ! H2o transmission over o3 band
1721 : real(r8), intent(out) :: abplnk1(14,pcols,pverp) ! non-nearest layer Plack factor
1722 : real(r8), intent(out) :: abplnk2(14,pcols,pverp) ! nearest layer factor
1723 :
1724 : !
1725 : !---------------------------Local variables-----------------------------
1726 : !
1727 : integer i ! Longitude index
1728 : integer k ! Level index]
1729 : integer k1 ! Level index
1730 : !
1731 : ! Local variables for H2O:
1732 : !
1733 : real(r8) h2oems(pcols,pverp) ! H2o emissivity
1734 : real(r8) tpathe ! Used to compute h2o emissivity
1735 : real(r8) dtx(pcols) ! Planck temperature minus 250 K
1736 : real(r8) dty(pcols) ! Path temperature minus 250 K
1737 : !
1738 : ! The 500-800 cm^-1 emission in emis(i,4) has been combined
1739 : ! into the 0-800 cm^-1 emission in emis(i,1)
1740 : !
1741 : real(r8) emis(pcols,2) ! H2O emissivity
1742 : !
1743 : !
1744 : !
1745 : real(r8) term7(pcols,2) ! Kl_inf(i) in eq(r8) of table A3a of R&D
1746 : real(r8) term8(pcols,2) ! Delta kl_inf(i) in eq(r8)
1747 : real(r8) tr1(pcols) ! Equation(6) in table A2 for 650-800
1748 : real(r8) tr2(pcols) ! Equation(6) in table A2 for 500-650
1749 : real(r8) tr3(pcols) ! Equation(4) in table A2 for 650-800
1750 : real(r8) tr4(pcols) ! Equation(4),table A2 of R&D for 500-650
1751 : real(r8) tr7(pcols) ! Equation (6) times eq(4) in table A2
1752 : ! of R&D for 650-800 cm-1 region
1753 : real(r8) tr8(pcols) ! Equation (6) times eq(4) in table A2
1754 : ! of R&D for 500-650 cm-1 region
1755 : real(r8) k21(pcols) ! Exponential coefficient used to calc
1756 : ! rot band transmissivity in the 650-800
1757 : ! cm-1 region (tr1)
1758 : real(r8) k22(pcols) ! Exponential coefficient used to calc
1759 : ! rot band transmissivity in the 500-650
1760 : ! cm-1 region (tr2)
1761 : real(r8) u(pcols) ! Pressure weighted H2O path length
1762 : real(r8) ub(nbands) ! Pressure weighted H2O path length with
1763 : ! Hulst-Curtis-Godson correction for
1764 : ! each band
1765 : real(r8) pnew ! Effective pressure for h2o linewidth
1766 : real(r8) pnewb(nbands) ! Effective pressure for h2o linewidth w/
1767 : ! Hulst-Curtis-Godson correction for
1768 : ! each band
1769 : real(r8) uc1(pcols) ! H2o continuum pathlength 500-800 cm-1
1770 : real(r8) fwk ! Equation(33) in R&D far wing correction
1771 : real(r8) troco2(pcols,pverp) ! H2o overlap factor for co2 absorption
1772 : real(r8) emplnk(14,pcols) ! emissivity Planck factor
1773 : real(r8) emstrc(pcols,pverp) ! total trace gas emissivity
1774 : !
1775 : ! Local variables for CO2:
1776 : !
1777 : real(r8) co2vmr(pcols) ! CO2 column mean vmr
1778 : real(r8) rmw ! ratio of molecular weights (air/co2)
1779 : real(r8) co2ems(pcols,pverp) ! Co2 emissivity
1780 : real(r8) co2plk(pcols) ! Used to compute co2 emissivity
1781 : real(r8) sum(pcols) ! Used to calculate path temperature
1782 : real(r8) t1i ! Co2 hot band temperature factor
1783 : real(r8) sqti ! Sqrt of temperature
1784 : real(r8) pi ! Pressure used in co2 mean line width
1785 : real(r8) et ! Co2 hot band factor
1786 : real(r8) et2 ! Co2 hot band factor
1787 : real(r8) et4 ! Co2 hot band factor
1788 : real(r8) omet ! Co2 stimulated emission term
1789 : real(r8) ex ! Part of co2 planck function
1790 : real(r8) f1co2 ! Co2 weak band factor
1791 : real(r8) f2co2 ! Co2 weak band factor
1792 : real(r8) f3co2 ! Co2 weak band factor
1793 : real(r8) t1co2 ! Overlap factor weak bands strong band
1794 : real(r8) sqwp ! Sqrt of co2 pathlength
1795 : real(r8) f1sqwp ! Main co2 band factor
1796 : real(r8) oneme ! Co2 stimulated emission term
1797 : real(r8) alphat ! Part of the co2 stimulated emiss term
1798 : real(r8) wco2 ! Consts used to define co2 pathlength
1799 : real(r8) posqt ! Effective pressure for co2 line width
1800 : real(r8) rbeta7 ! Inverse of co2 hot band line width par
1801 : real(r8) rbeta8 ! Inverse of co2 hot band line width par
1802 : real(r8) rbeta9 ! Inverse of co2 hot band line width par
1803 : real(r8) rbeta13 ! Inverse of co2 hot band line width par
1804 : real(r8) tpath ! Path temp used in co2 band model
1805 : real(r8) tmp1 ! Co2 band factor
1806 : real(r8) tmp2 ! Co2 band factor
1807 : real(r8) tmp3 ! Co2 band factor
1808 : real(r8) tlayr5 ! Temperature factor in co2 Planck func
1809 : real(r8) rsqti ! Reciprocal of sqrt of temperature
1810 : real(r8) exm1sq ! Part of co2 Planck function
1811 : real(r8) u7 ! Absorber amt for various co2 band systems
1812 : real(r8) u8 ! Absorber amt for various co2 band systems
1813 : real(r8) u9 ! Absorber amt for various co2 band systems
1814 : real(r8) u13 ! Absorber amt for various co2 band systems
1815 : real(r8) r250 ! Inverse 250K
1816 : real(r8) r300 ! Inverse 300K
1817 : real(r8) rsslp ! Inverse standard sea-level pressure
1818 : !
1819 : ! Local variables for O3:
1820 : !
1821 : real(r8) o3ems(pcols,pverp) ! Ozone emissivity
1822 : real(r8) dbvtt(pcols) ! Tmp drvtv of planck fctn for tplnke
1823 : real(r8) dbvt,fo3,t,ux,vx
1824 : real(r8) te ! Temperature factor
1825 : real(r8) u1 ! Path length factor
1826 : real(r8) u2 ! Path length factor
1827 : real(r8) phat ! Effecitive path length pressure
1828 : real(r8) tlocal ! Local planck function temperature
1829 : real(r8) tcrfac ! Scaled temperature factor
1830 : real(r8) beta ! Absorption funct factor voigt effect
1831 : real(r8) realnu ! Absorption function factor
1832 : real(r8) o3bndi ! Band absorption factor
1833 : !
1834 : ! Transmission terms for various spectral intervals:
1835 : !
1836 : real(r8) absbnd ! Proportional to co2 band absorptance
1837 : real(r8) tco2(pcols) ! co2 overlap factor
1838 : real(r8) th2o(pcols) ! h2o overlap factor
1839 : real(r8) to3(pcols) ! o3 overlap factor
1840 : !
1841 : ! Variables for new H2O parameterization
1842 : !
1843 : ! Notation:
1844 : ! U = integral (P/P_0 dW) eq. 15 in Ramanathan/Downey 1986
1845 : ! P = atmospheric pressure
1846 : ! P_0 = reference atmospheric pressure
1847 : ! W = precipitable water path
1848 : ! T_e = emission temperature
1849 : ! T_p = path temperature
1850 : ! RH = path relative humidity
1851 : !
1852 : real(r8) fe ! asymptotic value of emis. as U->infinity
1853 : real(r8) e_star ! normalized non-window emissivity
1854 : real(r8) l_star ! interpolated line transmission
1855 : real(r8) c_star ! interpolated continuum transmission
1856 :
1857 : real(r8) te1 ! emission temperature
1858 : real(r8) te2 ! te^2
1859 : real(r8) te3 ! te^3
1860 : real(r8) te4 ! te^4
1861 : real(r8) te5 ! te^5
1862 :
1863 : real(r8) log_u ! log base 10 of U
1864 : real(r8) log_uc ! log base 10 of H2O continuum path
1865 : real(r8) log_p ! log base 10 of P
1866 : real(r8) t_p ! T_p
1867 : real(r8) t_e ! T_e (offset by T_p)
1868 :
1869 : integer iu ! index for log10(U)
1870 : integer iu1 ! iu + 1
1871 : integer iuc ! index for log10(H2O continuum path)
1872 : integer iuc1 ! iuc + 1
1873 : integer ip ! index for log10(P)
1874 : integer ip1 ! ip + 1
1875 : integer itp ! index for T_p
1876 : integer itp1 ! itp + 1
1877 : integer ite ! index for T_e
1878 : integer ite1 ! ite + 1
1879 : integer irh ! index for RH
1880 : integer irh1 ! irh + 1
1881 :
1882 : real(r8) dvar ! normalized variation in T_p/T_e/P/U
1883 : real(r8) uvar ! U * diffusivity factor
1884 : real(r8) uscl ! factor for lineary scaling as U->0
1885 :
1886 : real(r8) wu ! weight for U
1887 : real(r8) wu1 ! 1 - wu
1888 : real(r8) wuc ! weight for H2O continuum path
1889 : real(r8) wuc1 ! 1 - wuc
1890 : real(r8) wp ! weight for P
1891 : real(r8) wp1 ! 1 - wp
1892 : real(r8) wtp ! weight for T_p
1893 : real(r8) wtp1 ! 1 - wtp
1894 : real(r8) wte ! weight for T_e
1895 : real(r8) wte1 ! 1 - wte
1896 : real(r8) wrh ! weight for RH
1897 : real(r8) wrh1 ! 1 - wrh
1898 :
1899 : real(r8) w_0_0_ ! weight for Tp/Te combination
1900 : real(r8) w_0_1_ ! weight for Tp/Te combination
1901 : real(r8) w_1_0_ ! weight for Tp/Te combination
1902 : real(r8) w_1_1_ ! weight for Tp/Te combination
1903 :
1904 : real(r8) w_0_00 ! weight for Tp/Te/RH combination
1905 : real(r8) w_0_01 ! weight for Tp/Te/RH combination
1906 : real(r8) w_0_10 ! weight for Tp/Te/RH combination
1907 : real(r8) w_0_11 ! weight for Tp/Te/RH combination
1908 : real(r8) w_1_00 ! weight for Tp/Te/RH combination
1909 : real(r8) w_1_01 ! weight for Tp/Te/RH combination
1910 : real(r8) w_1_10 ! weight for Tp/Te/RH combination
1911 : real(r8) w_1_11 ! weight for Tp/Te/RH combination
1912 :
1913 : real(r8) w00_00 ! weight for P/Tp/Te/RH combination
1914 : real(r8) w00_01 ! weight for P/Tp/Te/RH combination
1915 : real(r8) w00_10 ! weight for P/Tp/Te/RH combination
1916 : real(r8) w00_11 ! weight for P/Tp/Te/RH combination
1917 : real(r8) w01_00 ! weight for P/Tp/Te/RH combination
1918 : real(r8) w01_01 ! weight for P/Tp/Te/RH combination
1919 : real(r8) w01_10 ! weight for P/Tp/Te/RH combination
1920 : real(r8) w01_11 ! weight for P/Tp/Te/RH combination
1921 : real(r8) w10_00 ! weight for P/Tp/Te/RH combination
1922 : real(r8) w10_01 ! weight for P/Tp/Te/RH combination
1923 : real(r8) w10_10 ! weight for P/Tp/Te/RH combination
1924 : real(r8) w10_11 ! weight for P/Tp/Te/RH combination
1925 : real(r8) w11_00 ! weight for P/Tp/Te/RH combination
1926 : real(r8) w11_01 ! weight for P/Tp/Te/RH combination
1927 : real(r8) w11_10 ! weight for P/Tp/Te/RH combination
1928 : real(r8) w11_11 ! weight for P/Tp/Te/RH combination
1929 :
1930 : integer ib ! spectral interval:
1931 : ! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
1932 : ! 2 = 800-1200 cm^-1
1933 :
1934 : real(r8) pch2o ! H2O continuum path
1935 : real(r8) fch2o ! temp. factor for continuum
1936 : real(r8) uch2o ! U corresponding to H2O cont. path (window)
1937 :
1938 : real(r8) fdif ! secant(zenith angle) for diffusivity approx.
1939 :
1940 : real(r8) sslp_mks ! Sea-level pressure in MKS units
1941 : real(r8) esx ! saturation vapor pressure returned by qsat
1942 : real(r8) qsx ! saturation mixing ratio returned by qsat
1943 : real(r8) pnew_mks ! pnew in MKS units
1944 : real(r8) q_path ! effective specific humidity along path
1945 : real(r8) rh_path ! effective relative humidity along path
1946 :
1947 : !
1948 : !---------------------------Statement functions-------------------------
1949 : !
1950 : ! Derivative of planck function at 9.6 micro-meter wavelength, and
1951 : ! an absorption function factor:
1952 : !
1953 : !
1954 : dbvt(t)=(-2.8911366682e-4_r8+(2.3771251896e-6_r8+1.1305188929e-10_r8*t)*t)/ &
1955 : (1.0_r8+(-6.1364820707e-3_r8+1.5550319767e-5_r8*t)*t)
1956 : !
1957 : fo3(ux,vx)=ux/sqrt(4._r8+ux*(1._r8+vx))
1958 : !
1959 : !
1960 : !
1961 : !-----------------------------------------------------------------------
1962 : !
1963 : ! Initialize
1964 : !
1965 68112 : r250 = 1._r8/250._r8
1966 68112 : r300 = 1._r8/300._r8
1967 68112 : rsslp = 1._r8/sslp
1968 68112 : rmw = amd/amco2
1969 1137312 : do i=1,ncol
1970 1137312 : co2vmr(i) = co2mmr(i) * rmw
1971 : end do
1972 : !
1973 : ! Constants for computing U corresponding to H2O cont. path
1974 : !
1975 68112 : fdif = 1.66_r8
1976 68112 : sslp_mks = sslp / 10.0_r8
1977 : !
1978 : ! Planck function for co2
1979 : !
1980 1137312 : do i=1,ncol
1981 1069200 : ex = exp(960._r8/tplnke(i))
1982 1069200 : co2plk(i) = 5.e8_r8/((tplnke(i)**4)*(ex - 1._r8))
1983 1069200 : co2t(i,ntoplw) = tplnke(i)
1984 1137312 : sum(i) = co2t(i,ntoplw)*pnm(i,ntoplw)
1985 : end do
1986 68112 : k = ntoplw
1987 1839024 : do k1=pverp,ntoplw+1,-1
1988 1770912 : k = k + 1
1989 29638224 : do i=1,ncol
1990 27799200 : sum(i) = sum(i) + tlayr(i,k)*(pnm(i,k)-pnm(i,k-1))
1991 27799200 : ex = exp(960._r8/tlayr(i,k1))
1992 27799200 : tlayr5 = tlayr(i,k1)*tlayr4(i,k1)
1993 27799200 : co2eml(i,k1-1) = 1.2e11_r8*ex/(tlayr5*(ex - 1._r8)**2)
1994 29570112 : co2t(i,k) = sum(i)/pnm(i,k)
1995 : end do
1996 : end do
1997 : !
1998 : ! Initialize planck function derivative for O3
1999 : !
2000 1137312 : do i=1,ncol
2001 1137312 : dbvtt(i) = dbvt(tplnke(i))
2002 : end do
2003 : !
2004 : ! Calculate trace gas Planck functions
2005 : !
2006 : call trcplk(ncol , &
2007 : tint ,tlayr ,tplnke ,emplnk ,abplnk1 , &
2008 68112 : abplnk2 )
2009 :
2010 68112 : if ( ntoplw > 1 )then
2011 0 : emstot(:ncol,:ntoplw-1) = 0._r8
2012 : end if
2013 :
2014 : !
2015 : ! Interface loop
2016 : !
2017 1907136 : do k1=ntoplw,pverp
2018 : !
2019 : ! H2O emissivity
2020 : !
2021 : ! emis(i,1) 0 - 800 cm-1 h2o rotation band
2022 : ! emis(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
2023 : ! emis(i,2) 800 - 1200 cm-1 h2o window
2024 : !
2025 : ! Separation between rotation and vibration-rotation dropped, so
2026 : ! only 2 slots needed for H2O emissivity
2027 : !
2028 : ! emis(i,3) = 0.0
2029 : !
2030 : ! For the p type continuum
2031 : !
2032 30707424 : do i=1,ncol
2033 28868400 : u(i) = plh2o(i,k1)
2034 28868400 : pnew = u(i)/w(i,k1)
2035 28868400 : pnew_mks = pnew * sslp_mks
2036 : !
2037 : ! Apply scaling factor for 500-800 continuum
2038 : !
2039 : uc1(i) = (s2c(i,k1) + 1.7e-3_r8*plh2o(i,k1))*(1._r8 + 2._r8*s2c(i,k1))/ &
2040 28868400 : (1._r8 + 15._r8*s2c(i,k1))
2041 28868400 : pch2o = s2c(i,k1)
2042 : !
2043 : ! Changed effective path temperature to std. Curtis-Godson form
2044 : !
2045 28868400 : tpathe = tcg(i,k1)/w(i,k1)
2046 28868400 : t_p = min(max(tpathe, min_tp_h2o), max_tp_h2o)
2047 :
2048 28868400 : call qsat_water(t_p, pnew_mks, esx, qsx)
2049 :
2050 : !
2051 : ! Compute effective RH along path
2052 : !
2053 28868400 : q_path = w(i,k1) / pnm(i,k1) / rga
2054 : !
2055 : ! Calculate effective u, pnew for each band using
2056 : ! Hulst-Curtis-Godson approximation:
2057 : ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
2058 : ! 2nd edition, Oxford University Press, 1989.
2059 : ! Effective H2O path (w)
2060 : ! eq. 6.24, p. 228
2061 : ! Effective H2O path pressure (pnew = u/w):
2062 : ! eq. 6.29, p. 228
2063 : !
2064 28868400 : ub(1) = plh2ob(1,i,k1) / psi(t_p,1)
2065 28868400 : ub(2) = plh2ob(2,i,k1) / psi(t_p,2)
2066 :
2067 28868400 : pnewb(1) = ub(1) / wb(1,i,k1) * phi(t_p,1)
2068 28868400 : pnewb(2) = ub(2) / wb(2,i,k1) * phi(t_p,2)
2069 : !
2070 : !
2071 : !
2072 28868400 : dtx(i) = tplnke(i) - 250._r8
2073 28868400 : dty(i) = tpathe - 250._r8
2074 : !
2075 : ! Define variables for C/H/E (now C/LT/E) fit
2076 : !
2077 : ! emis(i,1) 0 - 800 cm-1 h2o rotation band
2078 : ! emis(i,1) 1200 - 2200 cm-1 h2o vibration-rotation band
2079 : ! emis(i,2) 800 - 1200 cm-1 h2o window
2080 : !
2081 : ! Separation between rotation and vibration-rotation dropped, so
2082 : ! only 2 slots needed for H2O emissivity
2083 : !
2084 : ! emis(i,3) = 0.0
2085 : !
2086 : ! Notation:
2087 : ! U = integral (P/P_0 dW)
2088 : ! P = atmospheric pressure
2089 : ! P_0 = reference atmospheric pressure
2090 : ! W = precipitable water path
2091 : ! T_e = emission temperature
2092 : ! T_p = path temperature
2093 : ! RH = path relative humidity
2094 : !
2095 : ! Terms for asymptotic value of emissivity
2096 : !
2097 28868400 : te1 = tplnke(i)
2098 28868400 : te2 = te1 * te1
2099 28868400 : te3 = te2 * te1
2100 28868400 : te4 = te3 * te1
2101 28868400 : te5 = te4 * te1
2102 : !
2103 : ! Band-independent indices for lines and continuum tables
2104 : !
2105 28868400 : dvar = (t_p - min_tp_h2o) / dtp_h2o
2106 28868400 : itp = min(max(int(aint(dvar,r8)) + 1, 1), n_tp - 1)
2107 28868400 : itp1 = itp + 1
2108 28868400 : wtp = dvar - floor(dvar)
2109 28868400 : wtp1 = 1.0_r8 - wtp
2110 :
2111 28868400 : t_e = min(max(tplnke(i) - t_p, min_te_h2o), max_te_h2o)
2112 28868400 : dvar = (t_e - min_te_h2o) / dte_h2o
2113 28868400 : ite = min(max(int(aint(dvar,r8)) + 1, 1), n_te - 1)
2114 28868400 : ite1 = ite + 1
2115 28868400 : wte = dvar - floor(dvar)
2116 28868400 : wte1 = 1.0_r8 - wte
2117 :
2118 28868400 : rh_path = min(max(q_path / qsx, min_rh_h2o), max_rh_h2o)
2119 28868400 : dvar = (rh_path - min_rh_h2o) / drh_h2o
2120 28868400 : irh = min(max(int(aint(dvar,r8)) + 1, 1), n_rh - 1)
2121 28868400 : irh1 = irh + 1
2122 28868400 : wrh = dvar - floor(dvar)
2123 28868400 : wrh1 = 1.0_r8 - wrh
2124 :
2125 28868400 : w_0_0_ = wtp * wte
2126 28868400 : w_0_1_ = wtp * wte1
2127 28868400 : w_1_0_ = wtp1 * wte
2128 28868400 : w_1_1_ = wtp1 * wte1
2129 :
2130 28868400 : w_0_00 = w_0_0_ * wrh
2131 28868400 : w_0_01 = w_0_0_ * wrh1
2132 28868400 : w_0_10 = w_0_1_ * wrh
2133 28868400 : w_0_11 = w_0_1_ * wrh1
2134 28868400 : w_1_00 = w_1_0_ * wrh
2135 28868400 : w_1_01 = w_1_0_ * wrh1
2136 28868400 : w_1_10 = w_1_1_ * wrh
2137 28868400 : w_1_11 = w_1_1_ * wrh1
2138 : !
2139 : ! H2O Continuum path for 0-800 and 1200-2200 cm^-1
2140 : !
2141 : ! Assume foreign continuum dominates total H2O continuum in these bands
2142 : ! per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776
2143 : ! Then the effective H2O path is just
2144 : ! U_c = integral[ f(P) dW ]
2145 : ! where
2146 : ! W = water-vapor mass and
2147 : ! f(P) = dependence of foreign continuum on pressure
2148 : ! = P / sslp
2149 : ! Then
2150 : ! U_c = U (the same effective H2O path as for lines)
2151 : !
2152 : !
2153 : ! Continuum terms for 800-1200 cm^-1
2154 : !
2155 : ! Assume self continuum dominates total H2O continuum for this band
2156 : ! per Clough et al, JGR, v. 97, no. D14 (Oct 20, 1992), p. 15776
2157 : ! Then the effective H2O self-continuum path is
2158 : ! U_c = integral[ h(e,T) dW ] (*eq. 1*)
2159 : ! where
2160 : ! W = water-vapor mass and
2161 : ! e = partial pressure of H2O along path
2162 : ! T = temperature along path
2163 : ! h(e,T) = dependence of foreign continuum on e,T
2164 : ! = e / sslp * f(T)
2165 : !
2166 : ! Replacing
2167 : ! e =~ q * P / epsilo
2168 : ! q = mixing ratio of H2O
2169 : ! epsilo = 0.622
2170 : !
2171 : ! and using the definition
2172 : ! U = integral [ (P / sslp) dW ]
2173 : ! = (P / sslp) W (homogeneous path)
2174 : !
2175 : ! the effective path length for the self continuum is
2176 : ! U_c = (q / epsilo) f(T) U (*eq. 2*)
2177 : !
2178 : ! Once values of T, U, and q have been calculated for the inhomogeneous
2179 : ! path, this sets U_c for the corresponding
2180 : ! homogeneous atmosphere. However, this need not equal the
2181 : ! value of U_c' defined by eq. 1 for the actual inhomogeneous atmosphere
2182 : ! under consideration.
2183 : !
2184 : ! Solution: hold T and q constant, solve for U' that gives U_c' by
2185 : ! inverting eq. (2):
2186 : !
2187 : ! U' = (U_c * epsilo) / (q * f(T))
2188 : !
2189 28868400 : fch2o = fh2oself(t_p)
2190 28868400 : uch2o = (pch2o * epsilo) / (q_path * fch2o)
2191 :
2192 : !
2193 : ! Band-dependent indices for non-window
2194 : !
2195 28868400 : ib = 1
2196 :
2197 28868400 : uvar = ub(ib) * fdif
2198 28868400 : log_u = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
2199 28868400 : dvar = (log_u - min_lu_h2o) / dlu_h2o
2200 28868400 : iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
2201 28868400 : iu1 = iu + 1
2202 28868400 : wu = dvar - floor(dvar)
2203 28868400 : wu1 = 1.0_r8 - wu
2204 :
2205 28868400 : log_p = min(log10(max(pnewb(ib), min_p_h2o)), max_lp_h2o)
2206 28868400 : dvar = (log_p - min_lp_h2o) / dlp_h2o
2207 28868400 : ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
2208 28868400 : ip1 = ip + 1
2209 28868400 : wp = dvar - floor(dvar)
2210 28868400 : wp1 = 1.0_r8 - wp
2211 :
2212 28868400 : w00_00 = wp * w_0_00
2213 28868400 : w00_01 = wp * w_0_01
2214 28868400 : w00_10 = wp * w_0_10
2215 28868400 : w00_11 = wp * w_0_11
2216 28868400 : w01_00 = wp * w_1_00
2217 28868400 : w01_01 = wp * w_1_01
2218 28868400 : w01_10 = wp * w_1_10
2219 28868400 : w01_11 = wp * w_1_11
2220 28868400 : w10_00 = wp1 * w_0_00
2221 28868400 : w10_01 = wp1 * w_0_01
2222 28868400 : w10_10 = wp1 * w_0_10
2223 28868400 : w10_11 = wp1 * w_0_11
2224 28868400 : w11_00 = wp1 * w_1_00
2225 28868400 : w11_01 = wp1 * w_1_01
2226 28868400 : w11_10 = wp1 * w_1_10
2227 28868400 : w11_11 = wp1 * w_1_11
2228 :
2229 : !
2230 : ! Asymptotic value of emissivity as U->infinity
2231 : !
2232 : fe = fet(1,ib) + &
2233 : fet(2,ib) * te1 + &
2234 : fet(3,ib) * te2 + &
2235 : fet(4,ib) * te3 + &
2236 : fet(5,ib) * te4 + &
2237 28868400 : fet(6,ib) * te5
2238 :
2239 : e_star = &
2240 144342000 : eh2onw(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
2241 28868400 : eh2onw(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
2242 28868400 : eh2onw(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
2243 : eh2onw(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
2244 28868400 : eh2onw(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
2245 : eh2onw(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
2246 : eh2onw(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
2247 : eh2onw(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
2248 28868400 : eh2onw(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
2249 : eh2onw(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
2250 : eh2onw(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
2251 : eh2onw(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
2252 : eh2onw(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
2253 : eh2onw(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
2254 : eh2onw(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
2255 : eh2onw(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
2256 28868400 : eh2onw(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
2257 : eh2onw(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
2258 : eh2onw(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
2259 : eh2onw(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
2260 : eh2onw(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
2261 : eh2onw(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
2262 : eh2onw(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
2263 : eh2onw(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
2264 : eh2onw(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
2265 : eh2onw(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
2266 : eh2onw(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
2267 : eh2onw(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
2268 : eh2onw(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
2269 : eh2onw(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
2270 : eh2onw(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
2271 288684000 : eh2onw(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
2272 : emis(i,ib) = min(max(fe * (1.0_r8 - (1.0_r8 - e_star) * &
2273 : aer_trn_ttl(i,k1,1,ib)), &
2274 28868400 : 0.0_r8), 1.0_r8)
2275 : !
2276 : ! Invoke linear limit for scaling wrt u below min_u_h2o
2277 : !
2278 28868400 : if (uvar < min_u_h2o) then
2279 3370444 : uscl = uvar / min_u_h2o
2280 3370444 : emis(i,ib) = emis(i,ib) * uscl
2281 : endif
2282 :
2283 :
2284 :
2285 : !
2286 : ! Band-dependent indices for window
2287 : !
2288 28868400 : ib = 2
2289 :
2290 28868400 : uvar = ub(ib) * fdif
2291 28868400 : log_u = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
2292 28868400 : dvar = (log_u - min_lu_h2o) / dlu_h2o
2293 28868400 : iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
2294 28868400 : iu1 = iu + 1
2295 28868400 : wu = dvar - floor(dvar)
2296 28868400 : wu1 = 1.0_r8 - wu
2297 :
2298 28868400 : log_p = min(log10(max(pnewb(ib), min_p_h2o)), max_lp_h2o)
2299 28868400 : dvar = (log_p - min_lp_h2o) / dlp_h2o
2300 28868400 : ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
2301 28868400 : ip1 = ip + 1
2302 28868400 : wp = dvar - floor(dvar)
2303 28868400 : wp1 = 1.0_r8 - wp
2304 :
2305 28868400 : w00_00 = wp * w_0_00
2306 28868400 : w00_01 = wp * w_0_01
2307 28868400 : w00_10 = wp * w_0_10
2308 28868400 : w00_11 = wp * w_0_11
2309 28868400 : w01_00 = wp * w_1_00
2310 28868400 : w01_01 = wp * w_1_01
2311 28868400 : w01_10 = wp * w_1_10
2312 28868400 : w01_11 = wp * w_1_11
2313 28868400 : w10_00 = wp1 * w_0_00
2314 28868400 : w10_01 = wp1 * w_0_01
2315 28868400 : w10_10 = wp1 * w_0_10
2316 28868400 : w10_11 = wp1 * w_0_11
2317 28868400 : w11_00 = wp1 * w_1_00
2318 28868400 : w11_01 = wp1 * w_1_01
2319 28868400 : w11_10 = wp1 * w_1_10
2320 28868400 : w11_11 = wp1 * w_1_11
2321 :
2322 28868400 : log_uc = min(log10(max(uch2o * fdif, min_u_h2o)), max_lu_h2o)
2323 28868400 : dvar = (log_uc - min_lu_h2o) / dlu_h2o
2324 28868400 : iuc = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
2325 28868400 : iuc1 = iuc + 1
2326 28868400 : wuc = dvar - floor(dvar)
2327 28868400 : wuc1 = 1.0_r8 - wuc
2328 : !
2329 : ! Asymptotic value of emissivity as U->infinity
2330 : !
2331 : fe = fet(1,ib) + &
2332 : fet(2,ib) * te1 + &
2333 : fet(3,ib) * te2 + &
2334 : fet(4,ib) * te3 + &
2335 : fet(5,ib) * te4 + &
2336 28868400 : fet(6,ib) * te5
2337 :
2338 : l_star = &
2339 57736800 : ln_eh2ow(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
2340 : ln_eh2ow(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
2341 : ln_eh2ow(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
2342 : ln_eh2ow(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
2343 28868400 : ln_eh2ow(ip , itp , iu1, ite , irh ) * w11_11 * wu + &
2344 : ln_eh2ow(ip , itp , iu1, ite , irh1) * w11_10 * wu + &
2345 : ln_eh2ow(ip , itp , iu1, ite1, irh ) * w11_01 * wu + &
2346 : ln_eh2ow(ip , itp , iu1, ite1, irh1) * w11_00 * wu + &
2347 : ln_eh2ow(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
2348 : ln_eh2ow(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
2349 : ln_eh2ow(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
2350 : ln_eh2ow(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
2351 : ln_eh2ow(ip , itp1, iu1, ite , irh ) * w10_11 * wu + &
2352 : ln_eh2ow(ip , itp1, iu1, ite , irh1) * w10_10 * wu + &
2353 : ln_eh2ow(ip , itp1, iu1, ite1, irh ) * w10_01 * wu + &
2354 : ln_eh2ow(ip , itp1, iu1, ite1, irh1) * w10_00 * wu + &
2355 28868400 : ln_eh2ow(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
2356 : ln_eh2ow(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
2357 : ln_eh2ow(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
2358 : ln_eh2ow(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
2359 : ln_eh2ow(ip1, itp , iu1, ite , irh ) * w01_11 * wu + &
2360 : ln_eh2ow(ip1, itp , iu1, ite , irh1) * w01_10 * wu + &
2361 : ln_eh2ow(ip1, itp , iu1, ite1, irh ) * w01_01 * wu + &
2362 : ln_eh2ow(ip1, itp , iu1, ite1, irh1) * w01_00 * wu + &
2363 : ln_eh2ow(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
2364 : ln_eh2ow(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
2365 : ln_eh2ow(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
2366 : ln_eh2ow(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
2367 : ln_eh2ow(ip1, itp1, iu1, ite , irh ) * w00_11 * wu + &
2368 : ln_eh2ow(ip1, itp1, iu1, ite , irh1) * w00_10 * wu + &
2369 : ln_eh2ow(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu + &
2370 115473600 : ln_eh2ow(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
2371 :
2372 : c_star = &
2373 28868400 : cn_eh2ow(ip , itp , iuc , ite , irh ) * w11_11 * wuc1 + &
2374 : cn_eh2ow(ip , itp , iuc , ite , irh1) * w11_10 * wuc1 + &
2375 : cn_eh2ow(ip , itp , iuc , ite1, irh ) * w11_01 * wuc1 + &
2376 : cn_eh2ow(ip , itp , iuc , ite1, irh1) * w11_00 * wuc1 + &
2377 28868400 : cn_eh2ow(ip , itp , iuc1, ite , irh ) * w11_11 * wuc + &
2378 : cn_eh2ow(ip , itp , iuc1, ite , irh1) * w11_10 * wuc + &
2379 : cn_eh2ow(ip , itp , iuc1, ite1, irh ) * w11_01 * wuc + &
2380 : cn_eh2ow(ip , itp , iuc1, ite1, irh1) * w11_00 * wuc + &
2381 : cn_eh2ow(ip , itp1, iuc , ite , irh ) * w10_11 * wuc1 + &
2382 : cn_eh2ow(ip , itp1, iuc , ite , irh1) * w10_10 * wuc1 + &
2383 : cn_eh2ow(ip , itp1, iuc , ite1, irh ) * w10_01 * wuc1 + &
2384 : cn_eh2ow(ip , itp1, iuc , ite1, irh1) * w10_00 * wuc1 + &
2385 : cn_eh2ow(ip , itp1, iuc1, ite , irh ) * w10_11 * wuc + &
2386 : cn_eh2ow(ip , itp1, iuc1, ite , irh1) * w10_10 * wuc + &
2387 : cn_eh2ow(ip , itp1, iuc1, ite1, irh ) * w10_01 * wuc + &
2388 : cn_eh2ow(ip , itp1, iuc1, ite1, irh1) * w10_00 * wuc + &
2389 : cn_eh2ow(ip1, itp , iuc , ite , irh ) * w01_11 * wuc1 + &
2390 : cn_eh2ow(ip1, itp , iuc , ite , irh1) * w01_10 * wuc1 + &
2391 : cn_eh2ow(ip1, itp , iuc , ite1, irh ) * w01_01 * wuc1 + &
2392 : cn_eh2ow(ip1, itp , iuc , ite1, irh1) * w01_00 * wuc1 + &
2393 : cn_eh2ow(ip1, itp , iuc1, ite , irh ) * w01_11 * wuc + &
2394 : cn_eh2ow(ip1, itp , iuc1, ite , irh1) * w01_10 * wuc + &
2395 : cn_eh2ow(ip1, itp , iuc1, ite1, irh ) * w01_01 * wuc + &
2396 : cn_eh2ow(ip1, itp , iuc1, ite1, irh1) * w01_00 * wuc + &
2397 : cn_eh2ow(ip1, itp1, iuc , ite , irh ) * w00_11 * wuc1 + &
2398 : cn_eh2ow(ip1, itp1, iuc , ite , irh1) * w00_10 * wuc1 + &
2399 : cn_eh2ow(ip1, itp1, iuc , ite1, irh ) * w00_01 * wuc1 + &
2400 : cn_eh2ow(ip1, itp1, iuc , ite1, irh1) * w00_00 * wuc1 + &
2401 : cn_eh2ow(ip1, itp1, iuc1, ite , irh ) * w00_11 * wuc + &
2402 : cn_eh2ow(ip1, itp1, iuc1, ite , irh1) * w00_10 * wuc + &
2403 : cn_eh2ow(ip1, itp1, iuc1, ite1, irh ) * w00_01 * wuc + &
2404 57736800 : cn_eh2ow(ip1, itp1, iuc1, ite1, irh1) * w00_00 * wuc
2405 : emis(i,ib) = min(max(fe * (1.0_r8 - l_star * c_star * &
2406 : aer_trn_ttl(i,k1,1,ib)), &
2407 28868400 : 0.0_r8), 1.0_r8)
2408 : !
2409 : ! Invoke linear limit for scaling wrt u below min_u_h2o
2410 : !
2411 28868400 : if (uvar < min_u_h2o) then
2412 3383658 : uscl = uvar / min_u_h2o
2413 3383658 : emis(i,ib) = emis(i,ib) * uscl
2414 : endif
2415 :
2416 :
2417 : !
2418 : ! Compute total emissivity for H2O
2419 : !
2420 59575824 : h2oems(i,k1) = emis(i,1)+emis(i,2)
2421 :
2422 : end do
2423 : !
2424 : !
2425 : !
2426 :
2427 30707424 : do i=1,ncol
2428 28868400 : term7(i,1) = coefj(1,1) + coefj(2,1)*dty(i)*(1._r8+c16*dty(i))
2429 28868400 : term8(i,1) = coefk(1,1) + coefk(2,1)*dty(i)*(1._r8+c17*dty(i))
2430 28868400 : term7(i,2) = coefj(1,2) + coefj(2,2)*dty(i)*(1._r8+c26*dty(i))
2431 30707424 : term8(i,2) = coefk(1,2) + coefk(2,2)*dty(i)*(1._r8+c27*dty(i))
2432 : end do
2433 30707424 : do i=1,ncol
2434 : !
2435 : ! 500 - 800 cm-1 rotation band overlap with co2
2436 : !
2437 28868400 : k21(i) = term7(i,1) + term8(i,1)/ &
2438 28868400 : (1._r8 + (c30 + c31*(dty(i)-10._r8)*(dty(i)-10._r8))*sqrt(u(i)))
2439 : k22(i) = term7(i,2) + term8(i,2)/ &
2440 28868400 : (1._r8 + (c28 + c29*(dty(i)-10._r8))*sqrt(u(i)))
2441 28868400 : fwk = fwcoef + fwc1/(1._r8+fwc2*u(i))
2442 28868400 : tr1(i) = exp(-(k21(i)*(sqrt(u(i)) + fc1*fwk*u(i))))
2443 28868400 : tr2(i) = exp(-(k22(i)*(sqrt(u(i)) + fc1*fwk*u(i))))
2444 28868400 : tr1(i)=tr1(i)*aer_trn_ttl(i,k1,1,idx_LW_0650_0800)
2445 : ! ! H2O line+aer trn 650--800 cm-1
2446 28868400 : tr2(i)=tr2(i)*aer_trn_ttl(i,k1,1,idx_LW_0500_0650)
2447 : ! ! H2O line+aer trn 500--650 cm-1
2448 28868400 : tr3(i) = exp(-((coefh(1,1) + coefh(2,1)*dtx(i))*uc1(i)))
2449 28868400 : tr4(i) = exp(-((coefh(1,2) + coefh(2,2)*dtx(i))*uc1(i)))
2450 28868400 : tr7(i) = tr1(i)*tr3(i)
2451 28868400 : tr8(i) = tr2(i)*tr4(i)
2452 28868400 : troco2(i,k1) = 0.65_r8*tr7(i) + 0.35_r8*tr8(i)
2453 30707424 : th2o(i) = tr8(i)
2454 : end do
2455 : !
2456 : ! CO2 emissivity for 15 micron band system
2457 : !
2458 30707424 : do i=1,ncol
2459 28868400 : t1i = exp(-480._r8/co2t(i,k1))
2460 28868400 : sqti = sqrt(co2t(i,k1))
2461 28868400 : rsqti = 1._r8/sqti
2462 28868400 : et = t1i
2463 28868400 : et2 = et*et
2464 28868400 : et4 = et2*et2
2465 28868400 : omet = 1._r8 - 1.5_r8*et2
2466 28868400 : f1co2 = 899.70_r8*omet*(1._r8 + 1.94774_r8*et + 4.73486_r8*et2)*rsqti
2467 28868400 : sqwp = sqrt(plco2(i,k1))
2468 28868400 : f1sqwp = f1co2*sqwp
2469 28868400 : t1co2 = 1._r8/(1._r8 + 245.18_r8*omet*sqwp*rsqti)
2470 28868400 : oneme = 1._r8 - et2
2471 28868400 : alphat = oneme**3*rsqti
2472 28868400 : wco2 = 2.5221_r8*co2vmr(i)*pnm(i,k1)*rga
2473 28868400 : u7 = 4.9411e4_r8*alphat*et2*wco2
2474 28868400 : u8 = 3.9744e4_r8*alphat*et4*wco2
2475 28868400 : u9 = 1.0447e5_r8*alphat*et4*et2*wco2
2476 28868400 : u13 = 2.8388e3_r8*alphat*et4*wco2
2477 : !
2478 28868400 : tpath = co2t(i,k1)
2479 28868400 : tlocal = tplnke(i)
2480 28868400 : tcrfac = sqrt((tlocal*r250)*(tpath*r300))
2481 28868400 : pi = pnm(i,k1)*rsslp + 2._r8*dpfco2*tcrfac
2482 28868400 : posqt = pi/(2._r8*sqti)
2483 28868400 : rbeta7 = 1._r8/( 5.3288_r8*posqt)
2484 28868400 : rbeta8 = 1._r8/ (10.6576_r8*posqt)
2485 28868400 : rbeta9 = rbeta7
2486 28868400 : rbeta13= rbeta9
2487 : f2co2 = (u7/sqrt(4._r8 + u7*(1._r8 + rbeta7))) + &
2488 : (u8/sqrt(4._r8 + u8*(1._r8 + rbeta8))) + &
2489 28868400 : (u9/sqrt(4._r8 + u9*(1._r8 + rbeta9)))
2490 28868400 : f3co2 = u13/sqrt(4._r8 + u13*(1._r8 + rbeta13))
2491 28868400 : tmp1 = log(1._r8 + f1sqwp)
2492 28868400 : tmp2 = log(1._r8 + f2co2)
2493 28868400 : tmp3 = log(1._r8 + f3co2)
2494 28868400 : absbnd = (tmp1 + 2._r8*t1co2*tmp2 + 2._r8*tmp3)*sqti
2495 28868400 : tco2(i)=1.0_r8/(1.0_r8+10.0_r8*(u7/sqrt(4._r8 + u7*(1._r8 + rbeta7))))
2496 28868400 : co2ems(i,k1) = troco2(i,k1)*absbnd*co2plk(i)
2497 28868400 : ex = exp(960._r8/tint(i,k1))
2498 28868400 : exm1sq = (ex - 1._r8)**2
2499 30707424 : co2em(i,k1) = 1.2e11_r8*ex/(tint(i,k1)*tint4(i,k1)*exm1sq)
2500 : end do
2501 : !
2502 : ! O3 emissivity
2503 : !
2504 30707424 : do i=1,ncol
2505 28868400 : h2otr(i,k1) = exp(-12._r8*s2c(i,k1))
2506 28868400 : h2otr(i,k1)=h2otr(i,k1)*aer_trn_ttl(i,k1,1,idx_LW_1000_1200)
2507 28868400 : te = (co2t(i,k1)/293._r8)**.7_r8
2508 28868400 : u1 = 18.29_r8*plos(i,k1)/te
2509 28868400 : u2 = .5649_r8*plos(i,k1)/te
2510 28868400 : phat = plos(i,k1)/plol(i,k1)
2511 28868400 : tlocal = tplnke(i)
2512 28868400 : tcrfac = sqrt(tlocal*r250)*te
2513 28868400 : beta = (1._r8/.3205_r8)*((1._r8/phat) + (dpfo3*tcrfac))
2514 28868400 : realnu = (1._r8/beta)*te
2515 28868400 : o3bndi = 74._r8*te*(tplnke(i)/375._r8)*log(1._r8 + fo3(u1,realnu) + fo3(u2,realnu))
2516 28868400 : o3ems(i,k1) = dbvtt(i)*h2otr(i,k1)*o3bndi
2517 30707424 : to3(i)=1.0_r8/(1._r8 + 0.1_r8*fo3(u1,realnu) + 0.1_r8*fo3(u2,realnu))
2518 : end do
2519 : !
2520 : ! Calculate trace gas emissivities
2521 : !
2522 : call trcems(ncol , &
2523 : k1 ,co2t ,pnm ,ucfc11 ,ucfc12 , &
2524 : un2o0 ,un2o1 ,bn2o0 ,bn2o1 ,uch4 , &
2525 : bch4 ,uco211 ,uco212 ,uco213 ,uco221 , &
2526 : uco222 ,uco223 ,uptype ,w ,s2c , &
2527 : u ,emplnk ,th2o ,tco2 ,to3 , &
2528 : emstrc , &
2529 1839024 : aer_trn_ttl)
2530 : !
2531 : ! Total emissivity:
2532 : !
2533 30775536 : do i=1,ncol
2534 57736800 : emstot(i,k1) = h2oems(i,k1) + co2ems(i,k1) + o3ems(i,k1) &
2535 88444224 : + emstrc(i,k1)
2536 : end do
2537 : end do ! End of interface loop
2538 :
2539 68112 : end subroutine radems
2540 :
2541 : !====================================================================================
2542 :
2543 749232 : subroutine radtpl(ncol , &
2544 : tnm ,lwupcgs ,qnm ,pnm ,plco2 ,plh2o , &
2545 : tplnka ,s2c ,tcg ,w ,tplnke , &
2546 : tint ,tint4 ,tlayr ,tlayr4 ,pmln , &
2547 : piln ,plh2ob ,wb ,co2mmr)
2548 : !--------------------------------------------------------------------
2549 : !
2550 : ! Purpose:
2551 : ! Compute temperatures and path lengths for longwave radiation
2552 : !
2553 : ! Method:
2554 : ! <Describe the algorithm(s) used in the routine.>
2555 : ! <Also include any applicable external references.>
2556 : !
2557 : ! Author: CCM1
2558 : !------------------------------Arguments-----------------------------
2559 : !
2560 : ! Input arguments
2561 : !
2562 : integer, intent(in) :: ncol ! number of atmospheric columns
2563 :
2564 : real(r8), intent(in) :: tnm(pcols,pver) ! Model level temperatures
2565 : real(r8), intent(in) :: lwupcgs(pcols) ! Surface longwave up flux
2566 : real(r8), intent(in) :: qnm(pcols,pver) ! Model level specific humidity
2567 : real(r8), intent(in) :: pnm(pcols,pverp) ! Pressure at model interfaces (dynes/cm2)
2568 : real(r8), intent(in) :: pmln(pcols,pver) ! Ln(pmidm1)
2569 : real(r8), intent(in) :: piln(pcols,pverp) ! Ln(pintm1)
2570 : real(r8), intent(in) :: co2mmr(pcols) ! co2 column mean mass mixing ratio
2571 : !
2572 : ! Output arguments
2573 : !
2574 : real(r8), intent(out) :: plco2(pcols,pverp) ! Pressure weighted co2 path
2575 : real(r8), intent(out) :: plh2o(pcols,pverp) ! Pressure weighted h2o path
2576 : real(r8), intent(out) :: tplnka(pcols,pverp) ! Level temperature from interface temperatures
2577 : real(r8), intent(out) :: s2c(pcols,pverp) ! H2o continuum path length
2578 : real(r8), intent(out) :: tcg(pcols,pverp) ! H2o-mass-wgted temp. (Curtis-Godson approx.)
2579 : real(r8), intent(out) :: w(pcols,pverp) ! H2o path length
2580 : real(r8), intent(out) :: tplnke(pcols) ! Equal to tplnka
2581 : real(r8), intent(out) :: tint(pcols,pverp) ! Layer interface temperature
2582 : real(r8), intent(out) :: tint4(pcols,pverp) ! Tint to the 4th power
2583 : real(r8), intent(out) :: tlayr(pcols,pverp) ! K-1 level temperature
2584 : real(r8), intent(out) :: tlayr4(pcols,pverp) ! Tlayr to the 4th power
2585 : real(r8), intent(out) :: plh2ob(nbands,pcols,pverp)! Pressure weighted h2o path with
2586 : ! Hulst-Curtis-Godson temp. factor
2587 : ! for H2O bands
2588 : real(r8), intent(out) :: wb(nbands,pcols,pverp) ! H2o path length with
2589 : ! Hulst-Curtis-Godson temp. factor
2590 : ! for H2O bands
2591 :
2592 : !
2593 : !---------------------------Local variables--------------------------
2594 : !
2595 : integer i ! Longitude index
2596 : integer k ! Level index
2597 : integer kp1 ! Level index + 1
2598 :
2599 : real(r8) repsil ! Inver ratio mol weight h2o to dry air
2600 : real(r8) dy ! Thickness of layer for tmp interp
2601 : real(r8) dpnm ! Pressure thickness of layer
2602 : real(r8) dpnmsq ! Prs squared difference across layer
2603 : real(r8) dw ! Increment in H2O path length
2604 : real(r8) dplh2o ! Increment in plh2o
2605 : real(r8) cpwpl ! Const in co2 mix ratio to path length conversn
2606 :
2607 : !--------------------------------------------------------------------
2608 : !
2609 749232 : repsil = 1._r8/epsilo
2610 : !
2611 : ! Compute co2 and h2o paths
2612 : !
2613 749232 : cpwpl = 0.5_r8/(gravit_cgs*p0)
2614 12510432 : do i=1,ncol
2615 11761200 : plh2o(i,ntoplw) = rgsslp*qnm(i,ntoplw)*pnm(i,ntoplw)*pnm(i,ntoplw)
2616 12510432 : plco2(i,ntoplw) = co2mmr(i)*cpwpl*pnm(i,ntoplw)*pnm(i,ntoplw)
2617 : end do
2618 20229264 : do k=ntoplw,pver
2619 326020464 : do i=1,ncol
2620 917373600 : plh2o(i,k+1) = plh2o(i,k) + rgsslp* &
2621 1223164800 : (pnm(i,k+1)**2 - pnm(i,k)**2)*qnm(i,k)
2622 325271232 : plco2(i,k+1) = co2mmr(i)*cpwpl*pnm(i,k+1)**2
2623 : end do
2624 : end do
2625 : !
2626 : ! Set the top and bottom intermediate level temperatures,
2627 : ! top level planck temperature and top layer temp**4.
2628 : !
2629 : ! Tint is lower interface temperature
2630 : ! (not available for bottom layer, so use ground temperature)
2631 : !
2632 12510432 : do i=1,ncol
2633 11761200 : tint4(i,pverp) = lwupcgs(i)/stebol_cgs
2634 11761200 : tint(i,pverp) = sqrt(sqrt(tint4(i,pverp)))
2635 11761200 : tplnka(i,ntoplw) = tnm(i,ntoplw)
2636 11761200 : tint(i,ntoplw) = tplnka(i,ntoplw)
2637 11761200 : tlayr4(i,ntoplw) = tplnka(i,ntoplw)**4
2638 12510432 : tint4(i,ntoplw) = tlayr4(i,ntoplw)
2639 : end do
2640 : !
2641 : ! Intermediate level temperatures are computed using temperature
2642 : ! at the full level below less dy*delta t,between the full level
2643 : !
2644 19480032 : do k=ntoplw+1,pver
2645 313510032 : do i=1,ncol
2646 294030000 : dy = (piln(i,k) - pmln(i,k))/(pmln(i,k-1) - pmln(i,k))
2647 294030000 : tint(i,k) = tnm(i,k) - dy*(tnm(i,k)-tnm(i,k-1))
2648 312760800 : tint4(i,k) = tint(i,k)**4
2649 : end do
2650 : end do
2651 : !
2652 : ! Now set the layer temp=full level temperatures and establish a
2653 : ! planck temperature for absorption (tplnka) which is the average
2654 : ! the intermediate level temperatures. Note that tplnka is not
2655 : ! equal to the full level temperatures.
2656 : !
2657 20229264 : do k=ntoplw+1,pverp
2658 326020464 : do i=1,ncol
2659 305791200 : tlayr(i,k) = tnm(i,k-1)
2660 305791200 : tlayr4(i,k) = tlayr(i,k)**4
2661 325271232 : tplnka(i,k) = .5_r8*(tint(i,k) + tint(i,k-1))
2662 : end do
2663 : end do
2664 : !
2665 : ! Calculate tplank for emissivity calculation.
2666 : ! Assume isothermal tplnke i.e. all levels=ttop.
2667 : !
2668 12510432 : do i=1,ncol
2669 11761200 : tplnke(i) = tplnka(i,ntoplw)
2670 12510432 : tlayr(i,ntoplw) = tint(i,ntoplw)
2671 : end do
2672 : !
2673 : ! Now compute h2o path fields:
2674 : !
2675 12510432 : do i=1,ncol
2676 : !
2677 : ! Changed effective path temperature to std. Curtis-Godson form
2678 : !
2679 11761200 : tcg(i,ntoplw) = rga*qnm(i,ntoplw)*pnm(i,ntoplw)*tnm(i,ntoplw)
2680 11761200 : w(i,ntoplw) = sslp * (plh2o(i,ntoplw)*2._r8) / pnm(i,ntoplw)
2681 : !
2682 : ! Hulst-Curtis-Godson scaling for H2O path
2683 : !
2684 11761200 : wb(1,i,ntoplw) = w(i,ntoplw) * phi(tnm(i,ntoplw),1)
2685 11761200 : wb(2,i,ntoplw) = w(i,ntoplw) * phi(tnm(i,ntoplw),2)
2686 : !
2687 : ! Hulst-Curtis-Godson scaling for effective pressure along H2O path
2688 : !
2689 11761200 : plh2ob(1,i,ntoplw) = plh2o(i,ntoplw) * psi(tnm(i,ntoplw),1)
2690 11761200 : plh2ob(2,i,ntoplw) = plh2o(i,ntoplw) * psi(tnm(i,ntoplw),2)
2691 :
2692 12510432 : s2c(i,ntoplw) = plh2o(i,ntoplw)*fh2oself(tnm(i,ntoplw))*qnm(i,ntoplw)*repsil
2693 : end do
2694 :
2695 20229264 : do k=ntoplw,pver
2696 326020464 : do i=1,ncol
2697 305791200 : dpnm = pnm(i,k+1) - pnm(i,k)
2698 305791200 : dpnmsq = pnm(i,k+1)**2 - pnm(i,k)**2
2699 305791200 : dw = rga*qnm(i,k)*dpnm
2700 305791200 : kp1 = k+1
2701 305791200 : w(i,kp1) = w(i,k) + dw
2702 : !
2703 : ! Hulst-Curtis-Godson scaling for H2O path
2704 : !
2705 305791200 : wb(1,i,kp1) = wb(1,i,k) + dw * phi(tnm(i,k),1)
2706 305791200 : wb(2,i,kp1) = wb(2,i,k) + dw * phi(tnm(i,k),2)
2707 : !
2708 : ! Hulst-Curtis-Godson scaling for effective pressure along H2O path
2709 : !
2710 305791200 : dplh2o = plh2o(i,kp1) - plh2o(i,k)
2711 :
2712 305791200 : plh2ob(1,i,kp1) = plh2ob(1,i,k) + dplh2o * psi(tnm(i,k),1)
2713 305791200 : plh2ob(2,i,kp1) = plh2ob(2,i,k) + dplh2o * psi(tnm(i,k),2)
2714 : !
2715 : ! Changed effective path temperature to std. Curtis-Godson form
2716 : !
2717 305791200 : tcg(i,kp1) = tcg(i,k) + dw*tnm(i,k)
2718 : s2c(i,kp1) = s2c(i,k) + rgsslp*dpnmsq*qnm(i,k)* &
2719 325271232 : fh2oself(tnm(i,k))*qnm(i,k)*repsil
2720 : end do
2721 : end do
2722 :
2723 749232 : end subroutine radtpl
2724 :
2725 : !====================================================================================
2726 :
2727 1536 : subroutine radae_init( &
2728 : gravx, epsilox, stebol, pstdx, mwdryx, &
2729 : mwco2x, mwo3x, absems_data)
2730 : !
2731 : ! Initialize radae module data
2732 : !
2733 : use pio, only: file_desc_t, var_desc_t, pio_inq_dimid, pio_inquire_dimension, &
2734 : pio_inquire_variable, pio_inq_varid, pio_get_var, pio_nowrite, &
2735 : pio_max_var_dims, pio_max_name, pio_closefile
2736 : use cam_pio_utils,only: cam_pio_openfile
2737 : use ioFileMod, only: getfil
2738 : !
2739 : ! Input variables
2740 : !
2741 : real(r8), intent(in) :: gravx ! Acceleration due to gravity (m/s**2)
2742 : real(r8), intent(in) :: epsilox ! Ratio of mol. wght of H2O to dry air
2743 : real(r8), intent(in) :: stebol ! Stefan-Boltzmann's constant (MKS)
2744 : real(r8), intent(in) :: pstdx ! Standard pressure (pascals)
2745 : real(r8), intent(in) :: mwdryx ! Molecular weight of dry air
2746 : real(r8), intent(in) :: mwco2x ! Molecular weight of carbon dioxide
2747 : real(r8), intent(in) :: mwo3x ! Molecular weight of ozone
2748 :
2749 : character(len=*) :: absems_data ! full pathname for time-invariant absorption dataset
2750 :
2751 : !
2752 : ! Variables for loading absorptivity/emissivity
2753 : !
2754 : type(file_desc_T) :: ncid_ae ! NetCDF file id for abs/ems file
2755 :
2756 : integer pdimid ! pressure dimension id
2757 : integer psize ! pressure dimension size
2758 :
2759 : integer tpdimid ! path temperature dimension id
2760 : integer tpsize ! path temperature size
2761 :
2762 : integer tedimid ! emission temperature dimension id
2763 : integer tesize ! emission temperature size
2764 :
2765 : integer udimid ! u (H2O path) dimension id
2766 : integer usize ! u (H2O path) dimension size
2767 :
2768 : integer rhdimid ! relative humidity dimension id
2769 : integer rhsize ! relative humidity dimension size
2770 :
2771 : type(var_desc_t) :: ah2onwid ! var. id for non-wndw abs.
2772 : type(var_desc_t) :: eh2onwid ! var. id for non-wndw ems.
2773 : type(var_desc_t) :: ah2owid ! var. id for wndw abs. (adjacent layers)
2774 : type(var_desc_t) :: cn_ah2owid ! var. id for continuum trans. for wndw abs.
2775 : type(var_desc_t) :: cn_eh2owid ! var. id for continuum trans. for wndw ems.
2776 : type(var_desc_t) :: ln_ah2owid ! var. id for line trans. for wndw abs.
2777 : type(var_desc_t) :: ln_eh2owid ! var. id for line trans. for wndw ems.
2778 :
2779 : character*(PIO_MAX_NAME) tmpname! dummy variable for var/dim names
2780 : character(len=256) locfn ! local filename
2781 : integer tmptype ! dummy variable for variable type
2782 : integer ndims ! number of dimensions
2783 : integer dims(PIO_MAX_VAR_DIMS) ! vector of dimension ids
2784 : integer natt ! number of attributes
2785 :
2786 : integer ierr ! ierr flag returned from pio (pio handles errors internally so it is not checked)
2787 : !
2788 : ! Constants to set
2789 : !
2790 1536 : gravit = gravx
2791 1536 : gravit_cgs = 100._r8*gravx
2792 1536 : rga = 1._r8/gravit_cgs
2793 1536 : epsilo = epsilox
2794 1536 : omeps = 1._r8 - epsilo
2795 1536 : sslp = 1.013250e6_r8
2796 1536 : stebol_cgs = 1.e3_r8*stebol
2797 1536 : rgsslp = 0.5_r8/(gravit_cgs*sslp)
2798 1536 : dpfo3 = 2.5e-3_r8
2799 1536 : dpfco2 = 5.0e-3_r8
2800 :
2801 1536 : p0 = pstdx*10.0_r8
2802 1536 : amd = mwdryx
2803 1536 : amco2 = mwco2x
2804 1536 : mwo3 = mwo3x
2805 : !
2806 : ! Coefficients for h2o emissivity and absorptivity for overlap of H2O
2807 : ! and trace gases.
2808 : !
2809 1536 : c16 = coefj(3,1)/coefj(2,1)
2810 1536 : c17 = coefk(3,1)/coefk(2,1)
2811 1536 : c26 = coefj(3,2)/coefj(2,2)
2812 1536 : c27 = coefk(3,2)/coefk(2,2)
2813 1536 : c28 = .5_r8
2814 1536 : c29 = .002053_r8
2815 1536 : c30 = .1_r8
2816 1536 : c31 = 3.0e-5_r8
2817 : !
2818 : ! Initialize further longwave constants referring to far wing
2819 : ! correction for overlap of H2O and trace gases; R&D refers to:
2820 : !
2821 : ! Ramanathan, V. and P.Downey, 1986: A Nonisothermal
2822 : ! Emissivity and Absorptivity Formulation for Water Vapor
2823 : ! Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
2824 : !
2825 1536 : fwcoef = .1_r8 ! See eq(33) R&D
2826 1536 : fwc1 = .30_r8 ! See eq(33) R&D
2827 1536 : fwc2 = 4.5_r8 ! See eq(33) and eq(34) in R&D
2828 1536 : fc1 = 2.6_r8 ! See eq(34) R&D
2829 :
2830 1536 : call getfil(absems_data, locfn)
2831 1536 : call cam_pio_openfile(ncid_ae, locfn, PIO_NOWRITE)
2832 :
2833 1536 : ierr = pio_inq_dimid(ncid_ae, 'p', pdimid)
2834 1536 : ierr = pio_inquire_dimension(ncid_ae, pdimid, len=psize)
2835 :
2836 1536 : ierr = pio_inq_dimid(ncid_ae, 'tp', tpdimid)
2837 1536 : ierr = pio_inquire_dimension(ncid_ae, tpdimid, len=tpsize)
2838 :
2839 1536 : ierr = pio_inq_dimid(ncid_ae, 'te', tedimid)
2840 1536 : ierr = pio_inquire_dimension(ncid_ae, tedimid, len=tesize)
2841 :
2842 1536 : ierr = pio_inq_dimid(ncid_ae, 'u', udimid)
2843 1536 : ierr = pio_inquire_dimension(ncid_ae, udimid, len=usize)
2844 :
2845 1536 : ierr = pio_inq_dimid(ncid_ae, 'rh', rhdimid)
2846 1536 : ierr = pio_inquire_dimension(ncid_ae, rhdimid, len=rhsize)
2847 :
2848 : if (psize /= n_p .or. &
2849 : tpsize /= n_tp .or. &
2850 : usize /= n_u .or. &
2851 1536 : tesize /= n_te .or. &
2852 : rhsize /= n_rh) then
2853 0 : call endrun ('RADAEINI: dimensions for abs/ems do not match internal def.')
2854 : endif
2855 :
2856 1536 : ierr = pio_inq_varid(ncid_ae, 'ah2onw', ah2onwid)
2857 1536 : ierr = pio_inq_varid(ncid_ae, 'eh2onw', eh2onwid)
2858 1536 : ierr = pio_inq_varid(ncid_ae, 'ah2ow', ah2owid)
2859 1536 : ierr = pio_inq_varid(ncid_ae, 'cn_ah2ow', cn_ah2owid)
2860 1536 : ierr = pio_inq_varid(ncid_ae, 'cn_eh2ow', cn_eh2owid)
2861 1536 : ierr = pio_inq_varid(ncid_ae, 'ln_ah2ow', ln_ah2owid)
2862 1536 : ierr = pio_inq_varid(ncid_ae, 'ln_eh2ow', ln_eh2owid)
2863 :
2864 1536 : ierr = pio_inquire_variable(ncid_ae, ah2onwid, tmpname, tmptype, ndims, dims, natt)
2865 : if (ndims /= 5 .or. &
2866 : dims(1) /= pdimid .or. &
2867 : dims(2) /= tpdimid .or. &
2868 : dims(3) /= udimid .or. &
2869 1536 : dims(4) /= tedimid .or. &
2870 : dims(5) /= rhdimid) then
2871 0 : call endrun ('RADAEINI: non-wndw abs. in file /= internal def.')
2872 : endif
2873 1536 : ierr = pio_inquire_variable(ncid_ae, eh2onwid, tmpname, tmptype, ndims, dims, natt)
2874 : if (ndims /= 5 .or. &
2875 : dims(1) /= pdimid .or. &
2876 : dims(2) /= tpdimid .or. &
2877 : dims(3) /= udimid .or. &
2878 1536 : dims(4) /= tedimid .or. &
2879 : dims(5) /= rhdimid) then
2880 0 : call endrun ('RADAEINI: non-wndw ems. in file /= internal def.')
2881 : endif
2882 1536 : ierr = pio_inquire_variable(ncid_ae, ah2owid, tmpname, tmptype, ndims, dims, natt)
2883 : if (ndims /= 5 .or. &
2884 : dims(1) /= pdimid .or. &
2885 : dims(2) /= tpdimid .or. &
2886 : dims(3) /= udimid .or. &
2887 1536 : dims(4) /= tedimid .or. &
2888 : dims(5) /= rhdimid) then
2889 0 : call endrun ('RADAEINI: window abs. in file /= internal def.')
2890 : endif
2891 1536 : ierr = pio_inquire_variable(ncid_ae, cn_ah2owid, tmpname, tmptype, ndims, dims, natt)
2892 : if (ndims /= 5 .or. &
2893 : dims(1) /= pdimid .or. &
2894 : dims(2) /= tpdimid .or. &
2895 : dims(3) /= udimid .or. &
2896 1536 : dims(4) /= tedimid .or. &
2897 : dims(5) /= rhdimid) then
2898 0 : call endrun ('RADAEINI: cont. trans for abs. in file /= internal def.')
2899 : endif
2900 1536 : ierr = pio_inquire_variable(ncid_ae, cn_eh2owid, tmpname, tmptype, ndims, dims, natt)
2901 : if (ndims /= 5 .or. &
2902 : dims(1) /= pdimid .or. &
2903 : dims(2) /= tpdimid .or. &
2904 : dims(3) /= udimid .or. &
2905 1536 : dims(4) /= tedimid .or. &
2906 : dims(5) /= rhdimid) then
2907 0 : call endrun ('RADAEINI: cont. trans. for ems. in file /= internal def.')
2908 : endif
2909 1536 : ierr = pio_inquire_variable(ncid_ae, ln_ah2owid, tmpname, tmptype, ndims, dims, natt)
2910 : if (ndims /= 5 .or. &
2911 : dims(1) /= pdimid .or. &
2912 : dims(2) /= tpdimid .or. &
2913 : dims(3) /= udimid .or. &
2914 1536 : dims(4) /= tedimid .or. &
2915 : dims(5) /= rhdimid) then
2916 0 : call endrun ('RADAEINI: line trans for abs. in file /= internal def.')
2917 : endif
2918 1536 : ierr = pio_inquire_variable(ncid_ae, ln_eh2owid, tmpname, tmptype, ndims, dims, natt)
2919 : if (ndims /= 5 .or. &
2920 : dims(1) /= pdimid .or. &
2921 : dims(2) /= tpdimid .or. &
2922 : dims(3) /= udimid .or. &
2923 1536 : dims(4) /= tedimid .or. &
2924 : dims(5) /= rhdimid) then
2925 0 : call endrun ('RADAEINI: line trans. for ems. in file /= internal def.')
2926 : endif
2927 :
2928 1536 : ierr = pio_get_var (ncid_ae, ah2onwid, ah2onw)
2929 1536 : ierr = pio_get_var (ncid_ae, eh2onwid, eh2onw)
2930 1536 : ierr = pio_get_var (ncid_ae, ah2owid, ah2ow)
2931 1536 : ierr = pio_get_var (ncid_ae, cn_ah2owid, cn_ah2ow)
2932 1536 : ierr = pio_get_var (ncid_ae, cn_eh2owid, cn_eh2ow)
2933 1536 : ierr = pio_get_var (ncid_ae, ln_ah2owid, ln_ah2ow)
2934 1536 : ierr = pio_get_var (ncid_ae, ln_eh2owid, ln_eh2ow)
2935 :
2936 1536 : call pio_closefile(ncid_ae)
2937 :
2938 : ! check whether arrays have already been allocated before calling
2939 : ! initialize_radbuffer. This will be the case for restart and branch
2940 : ! runs since those arrays are restored from the restart file before
2941 : ! this routine is called.
2942 1536 : if (.not. allocated(abstot_3d)) call initialize_radbuffer()
2943 :
2944 3072 : end subroutine radae_init
2945 :
2946 : !====================================================================================
2947 :
2948 1536 : subroutine initialize_radbuffer
2949 : !
2950 : ! Initialize radiation buffer data
2951 : !
2952 :
2953 1536 : use ref_pres, only : pref_mid
2954 : use phys_control, only : phys_getopts
2955 :
2956 : character(len=16) :: radiation_scheme
2957 : integer :: k
2958 :
2959 : ! If the top model level is above ~90 km (0.1 Pa), set the top level to compute
2960 : ! longwave cooling to about 80 km (1 Pa)
2961 1536 : if (pref_mid(1) .lt. 0.1_r8) then
2962 0 : do k = 1, plev
2963 0 : if (pref_mid(k) .lt. 1._r8) ntoplw = k
2964 : end do
2965 : else
2966 1536 : ntoplw = 1
2967 : end if
2968 1536 : if (masterproc) then
2969 2 : write(iulog,*) 'INITIALIZE_RADBUFFER: ntoplw =',ntoplw, ' pressure:',pref_mid(ntoplw)
2970 : endif
2971 :
2972 1536 : call phys_getopts(radiation_scheme_out=radiation_scheme)
2973 :
2974 1536 : if(radiation_scheme.eq.'camrt') then
2975 7680 : allocate (abstot_3d(pcols,ntoplw:pverp,ntoplw:pverp,begchunk:endchunk))
2976 4608 : allocate (absnxt_3d(pcols,pver,4,begchunk:endchunk))
2977 4608 : allocate (emstot_3d(pcols,pverp,begchunk:endchunk))
2978 1536 : abstot_3d(:,:,:,:) = posinf
2979 1536 : absnxt_3d(:,:,:,:) = posinf
2980 1536 : emstot_3d(:,:,:) = posinf
2981 : end if
2982 1536 : return
2983 1536 : end subroutine initialize_radbuffer
2984 :
2985 : !====================================================================================
2986 :
2987 68112 : subroutine radoz2(ncol, o3, pint, plol, plos)
2988 : !-----------------------------------------------------------------------
2989 : !
2990 : ! Purpose:
2991 : ! Computes the path length integrals to the model interfaces given the
2992 : ! ozone volume mixing ratio
2993 : !
2994 : ! Method:
2995 : ! <Describe the algorithm(s) used in the routine.>
2996 : ! <Also include any applicable external references.>
2997 : !
2998 : ! Author: CCM1, CMS Contact J. Kiehl
2999 : !
3000 : !------------------------------Input arguments--------------------------
3001 : !
3002 : integer, intent(in) :: ncol ! number of atmospheric columns
3003 :
3004 : real(r8), intent(in) :: o3(pcols,pver) ! ozone mass mixing ratio
3005 : real(r8), intent(in) :: pint(pcols,pverp) ! Model interface pressures
3006 : !
3007 : !----------------------------Output arguments---------------------------
3008 : !
3009 : real(r8), intent(out) :: plol(pcols,pverp) ! Ozone prs weighted path length (cm)
3010 : real(r8), intent(out) :: plos(pcols,pverp) ! Ozone path length (cm)
3011 : !
3012 : !---------------------------Local workspace-----------------------------
3013 : !
3014 : integer i ! longitude index
3015 : integer k ! level index
3016 :
3017 : real(r8) :: v0 ! Volume of a gas at stp (m**3/kmol)
3018 : real(r8) :: p0 ! Standard pressure (pascals)
3019 : real(r8) :: cplos ! constant for ozone path length integral
3020 : real(r8) :: cplol ! constant for ozone path length integral
3021 :
3022 : !
3023 : !-----------------------------------------------------------------------
3024 : !*******************************************************************
3025 : ! These hardwired constants need to be replaced with common values.
3026 : ! They are here for testing infrastructure changes that should not
3027 : ! change answers.
3028 : ! Constants for ozone path integrals (multiplication by 100 for unit
3029 : ! conversion to cgs from mks):
3030 : !
3031 68112 : v0 = 22.4136_r8 ! Volume of a gas at stp (m**3/kmol)
3032 68112 : p0 = 0.1_r8*sslp ! Standard pressure (pascals)
3033 68112 : cplos = v0/(mwo3*gravit) *100.0_r8
3034 68112 : cplol = v0/(mwo3*gravit*p0)*0.5_r8*100.0_r8
3035 : !*******************************************************************
3036 : !
3037 : ! Evaluate the ozone path length integrals to interfaces;
3038 : ! factors of .1 and .01 to convert pressures from cgs to mks:
3039 : !
3040 1137312 : do i=1,ncol
3041 1069200 : plos(i,ntoplw) = 0.1_r8 *cplos*o3(i,ntoplw)*pint(i,ntoplw)
3042 1137312 : plol(i,ntoplw) = 0.01_r8*cplol*o3(i,ntoplw)*pint(i,ntoplw)*pint(i,ntoplw)
3043 : end do
3044 1839024 : do k=ntoplw+1,pverp
3045 29638224 : do i=1,ncol
3046 27799200 : plos(i,k) = plos(i,k-1) + 0.1_r8*cplos*o3(i,k-1)*(pint(i,k) - pint(i,k-1))
3047 : plol(i,k) = plol(i,k-1) + 0.01_r8*cplol*o3(i,k-1)* &
3048 29570112 : (pint(i,k)*pint(i,k) - pint(i,k-1)*pint(i,k-1))
3049 : end do
3050 : end do
3051 :
3052 1536 : end subroutine radoz2
3053 :
3054 : !====================================================================================
3055 :
3056 68112 : subroutine trcpth(ncol , &
3057 : tnm ,pnm ,cfc11 ,cfc12 ,n2o , &
3058 : ch4 ,qnm ,ucfc11 ,ucfc12 ,un2o0 , &
3059 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
3060 : uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
3061 : bch4 ,uptype ,co2mmr)
3062 : !-----------------------------------------------------------------------
3063 : !
3064 : ! Purpose:
3065 : ! Calculate path lengths and pressure factors for CH4, N2O, CFC11
3066 : ! and CFC12.
3067 : !
3068 : ! Method:
3069 : ! See CCM3 description for details
3070 : !
3071 : ! Author: J. Kiehl
3072 : !
3073 : !-----------------------------------------------------------------------
3074 : !
3075 : ! Input arguments
3076 : !
3077 : integer, intent(in) :: ncol ! number of atmospheric columns
3078 :
3079 : real(r8), intent(in) :: tnm(pcols,pver) ! Model level temperatures
3080 : real(r8), intent(in) :: pnm(pcols,pverp) ! Pres. at model interfaces (dynes/cm2)
3081 : real(r8), intent(in) :: qnm(pcols,pver) ! h2o specific humidity
3082 : real(r8), intent(in) :: cfc11(pcols,pver) ! CFC11 mass mixing ratio
3083 : !
3084 : real(r8), intent(in) :: cfc12(pcols,pver) ! CFC12 mass mixing ratio
3085 : real(r8), intent(in) :: n2o(pcols,pver) ! N2O mass mixing ratio
3086 : real(r8), intent(in) :: ch4(pcols,pver) ! CH4 mass mixing ratio
3087 : real(r8), intent(in) :: co2mmr(pcols) ! co2 column mean mass mixing ratio
3088 :
3089 : !
3090 : ! Output arguments
3091 : !
3092 : real(r8), intent(out) :: ucfc11(pcols,pverp) ! CFC11 path length
3093 : real(r8), intent(out) :: ucfc12(pcols,pverp) ! CFC12 path length
3094 : real(r8), intent(out) :: un2o0(pcols,pverp) ! N2O path length
3095 : real(r8), intent(out) :: un2o1(pcols,pverp) ! N2O path length (hot band)
3096 : real(r8), intent(out) :: uch4(pcols,pverp) ! CH4 path length
3097 : !
3098 : real(r8), intent(out) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
3099 : real(r8), intent(out) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
3100 : real(r8), intent(out) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
3101 : real(r8), intent(out) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
3102 : real(r8), intent(out) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
3103 : !
3104 : real(r8), intent(out) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
3105 : real(r8), intent(out) :: bn2o0(pcols,pverp) ! pressure factor for n2o
3106 : real(r8), intent(out) :: bn2o1(pcols,pverp) ! pressure factor for n2o
3107 : real(r8), intent(out) :: bch4(pcols,pverp) ! pressure factor for ch4
3108 : real(r8), intent(out) :: uptype(pcols,pverp) ! p-type continuum path length
3109 :
3110 : !
3111 : !---------------------------Local variables-----------------------------
3112 : !
3113 : integer i ! Longitude index
3114 : integer k ! Level index
3115 : !
3116 : real(r8) co2fac(pcols,1) ! co2 factor
3117 : real(r8) alpha1(pcols) ! stimulated emission term
3118 : real(r8) alpha2(pcols) ! stimulated emission term
3119 : real(r8) rt(pcols) ! reciprocal of local temperature
3120 : real(r8) rsqrt(pcols) ! reciprocal of sqrt of temp
3121 : !
3122 : real(r8) pbar(pcols) ! mean pressure
3123 : real(r8) dpnm(pcols) ! difference in pressure
3124 : real(r8) diff ! diffusivity factor
3125 : !
3126 : !--------------------------Data Statements------------------------------
3127 : !
3128 : data diff /1.66_r8/
3129 : !
3130 : !-----------------------------------------------------------------------
3131 : !
3132 : ! Calculate path lengths for the trace gases at model top
3133 : !
3134 :
3135 1137312 : do i = 1,ncol
3136 1069200 : ucfc11(i,ntoplw) = 1.8_r8 * cfc11(i,ntoplw) * pnm(i,ntoplw) * rga
3137 1069200 : ucfc12(i,ntoplw) = 1.8_r8 * cfc12(i,ntoplw) * pnm(i,ntoplw) * rga
3138 1069200 : un2o0(i,ntoplw) = diff * 1.02346e5_r8 * n2o(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
3139 1069200 : un2o1(i,ntoplw) = diff * 2.01909_r8 * un2o0(i,ntoplw) * exp(-847.36_r8/tnm(i,ntoplw))
3140 1069200 : uch4(i,ntoplw) = diff * 8.60957e4_r8 * ch4(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
3141 1069200 : co2fac(i,1) = diff * co2mmr(i) * pnm(i,ntoplw) * rga
3142 1069200 : alpha1(i) = (1.0_r8 - exp(-1540.0_r8/tnm(i,ntoplw)))**3.0_r8/sqrt(tnm(i,ntoplw))
3143 1069200 : alpha2(i) = (1.0_r8 - exp(-1360.0_r8/tnm(i,ntoplw)))**3.0_r8/sqrt(tnm(i,ntoplw))
3144 1069200 : uco211(i,ntoplw) = 3.42217e3_r8 * co2fac(i,1) * alpha1(i) * exp(-1849.7_r8/tnm(i,ntoplw))
3145 1069200 : uco212(i,ntoplw) = 6.02454e3_r8 * co2fac(i,1) * alpha1(i) * exp(-2782.1_r8/tnm(i,ntoplw))
3146 1069200 : uco213(i,ntoplw) = 5.53143e3_r8 * co2fac(i,1) * alpha1(i) * exp(-3723.2_r8/tnm(i,ntoplw))
3147 1069200 : uco221(i,ntoplw) = 3.88984e3_r8 * co2fac(i,1) * alpha2(i) * exp(-1997.6_r8/tnm(i,ntoplw))
3148 1069200 : uco222(i,ntoplw) = 3.67108e3_r8 * co2fac(i,1) * alpha2(i) * exp(-3843.8_r8/tnm(i,ntoplw))
3149 1069200 : uco223(i,ntoplw) = 6.50642e3_r8 * co2fac(i,1) * alpha2(i) * exp(-2989.7_r8/tnm(i,ntoplw))
3150 : bn2o0(i,ntoplw) = diff * 19.399_r8 * pnm(i,ntoplw)**2.0_r8 * n2o(i,ntoplw) * &
3151 1069200 : 1.02346e5_r8 * rga / (sslp*tnm(i,ntoplw))
3152 1069200 : bn2o1(i,ntoplw) = bn2o0(i,ntoplw) * exp(-847.36_r8/tnm(i,ntoplw)) * 2.06646e5_r8
3153 : bch4(i,ntoplw) = diff * 2.94449_r8 * ch4(i,ntoplw) * pnm(i,ntoplw)**2.0_r8 * rga * &
3154 1069200 : 8.60957e4_r8 / (sslp*tnm(i,ntoplw))
3155 : uptype(i,ntoplw) = diff * qnm(i,ntoplw) * pnm(i,ntoplw)**2.0_r8 * &
3156 1137312 : exp(1800.0_r8*(1.0_r8/tnm(i,ntoplw) - 1.0_r8/296.0_r8)) * rga / sslp
3157 : end do
3158 : !
3159 : ! Calculate trace gas path lengths through model atmosphere
3160 : !
3161 1839024 : do k = ntoplw,pver
3162 29638224 : do i = 1,ncol
3163 27799200 : rt(i) = 1._r8/tnm(i,k)
3164 27799200 : rsqrt(i) = sqrt(rt(i))
3165 27799200 : pbar(i) = 0.5_r8 * (pnm(i,k+1) + pnm(i,k)) / sslp
3166 27799200 : dpnm(i) = (pnm(i,k+1) - pnm(i,k)) * rga
3167 27799200 : alpha1(i) = diff * rsqrt(i) * (1.0_r8 - exp(-1540.0_r8/tnm(i,k)))**3.0_r8
3168 27799200 : alpha2(i) = diff * rsqrt(i) * (1.0_r8 - exp(-1360.0_r8/tnm(i,k)))**3.0_r8
3169 27799200 : ucfc11(i,k+1) = ucfc11(i,k) + 1.8_r8 * cfc11(i,k) * dpnm(i)
3170 27799200 : ucfc12(i,k+1) = ucfc12(i,k) + 1.8_r8 * cfc12(i,k) * dpnm(i)
3171 27799200 : un2o0(i,k+1) = un2o0(i,k) + diff * 1.02346e5_r8 * n2o(i,k) * rsqrt(i) * dpnm(i)
3172 : un2o1(i,k+1) = un2o1(i,k) + diff * 2.06646e5_r8 * n2o(i,k) * &
3173 27799200 : rsqrt(i) * exp(-847.36_r8/tnm(i,k)) * dpnm(i)
3174 27799200 : uch4(i,k+1) = uch4(i,k) + diff * 8.60957e4_r8 * ch4(i,k) * rsqrt(i) * dpnm(i)
3175 : uco211(i,k+1) = uco211(i,k) + 1.15_r8*3.42217e3_r8 * alpha1(i) * &
3176 27799200 : co2mmr(i) * exp(-1849.7_r8/tnm(i,k)) * dpnm(i)
3177 : uco212(i,k+1) = uco212(i,k) + 1.15_r8*6.02454e3_r8 * alpha1(i) * &
3178 27799200 : co2mmr(i) * exp(-2782.1_r8/tnm(i,k)) * dpnm(i)
3179 : uco213(i,k+1) = uco213(i,k) + 1.15_r8*5.53143e3_r8 * alpha1(i) * &
3180 27799200 : co2mmr(i) * exp(-3723.2_r8/tnm(i,k)) * dpnm(i)
3181 : uco221(i,k+1) = uco221(i,k) + 1.15_r8*3.88984e3_r8 * alpha2(i) * &
3182 27799200 : co2mmr(i) * exp(-1997.6_r8/tnm(i,k)) * dpnm(i)
3183 : uco222(i,k+1) = uco222(i,k) + 1.15_r8*3.67108e3_r8 * alpha2(i) * &
3184 27799200 : co2mmr(i) * exp(-3843.8_r8/tnm(i,k)) * dpnm(i)
3185 : uco223(i,k+1) = uco223(i,k) + 1.15_r8*6.50642e3_r8 * alpha2(i) * &
3186 27799200 : co2mmr(i) * exp(-2989.7_r8/tnm(i,k)) * dpnm(i)
3187 : bn2o0(i,k+1) = bn2o0(i,k) + diff * 19.399_r8 * pbar(i) * rt(i) &
3188 27799200 : * 1.02346e5_r8 * n2o(i,k) * dpnm(i)
3189 : bn2o1(i,k+1) = bn2o1(i,k) + diff * 19.399_r8 * pbar(i) * rt(i) &
3190 27799200 : * 2.06646e5_r8 * exp(-847.36_r8/tnm(i,k)) * n2o(i,k)*dpnm(i)
3191 : bch4(i,k+1) = bch4(i,k) + diff * 2.94449_r8 * rt(i) * pbar(i) &
3192 27799200 : * 8.60957e4_r8 * ch4(i,k) * dpnm(i)
3193 : uptype(i,k+1) = uptype(i,k) + diff *qnm(i,k) * &
3194 29570112 : exp(1800.0_r8*(1.0_r8/tnm(i,k) - 1.0_r8/296.0_r8)) * pbar(i) * dpnm(i)
3195 : end do
3196 : end do
3197 : !
3198 68112 : return
3199 : end subroutine trcpth
3200 :
3201 :
3202 :
3203 : !====================================================================================
3204 : ! Private Interfaces
3205 : !====================================================================================
3206 :
3207 1096999200 : function fh2oself( temp )
3208 : !
3209 : ! Short function for H2O self-continuum temperature factor in
3210 : ! calculation of effective H2O self-continuum path length
3211 : !
3212 : ! H2O Continuum: CKD 2.4
3213 : ! Code for continuum: GENLN3
3214 : ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
3215 : ! Transmittance and Radiance Model, Version 3.0 Description
3216 : ! and Users Guide, NCAR/TN-367+STR, 147 pp.
3217 : !
3218 : ! In GENLN, the temperature scaling of the self-continuum is handled
3219 : ! by exponential interpolation/extrapolation from observations at
3220 : ! 260K and 296K by:
3221 : !
3222 : ! TFAC = (T(IPATH) - 296.0)/(260.0 - 296.0)
3223 : ! CSFFT = CSFF296*(CSFF260/CSFF296)**TFAC
3224 : !
3225 : ! For 800-1200 cm^-1, (CSFF260/CSFF296) ranges from ~2.1 to ~1.9
3226 : ! with increasing wavenumber. The ratio <CSFF260>/<CSFF296>,
3227 : ! where <> indicates average over wavenumber, is ~2.07
3228 : !
3229 : ! fh2oself is (<CSFF260>/<CSFF296>)**TFAC
3230 : !
3231 : real(r8),intent(in) :: temp ! path temperature
3232 : real(r8) fh2oself ! mean ratio of self-continuum at temp and 296K
3233 :
3234 1096999200 : fh2oself = 2.0727484_r8**((296.0_r8 - temp) / 36.0_r8)
3235 1096999200 : end function fh2oself
3236 :
3237 : !====================================================================================
3238 :
3239 2193998400 : function phi(tpx,iband)
3240 : !
3241 : ! History: First version for Hitran 1996 (C/H/E)
3242 : ! Current version for Hitran 2000 (C/LT/E)
3243 : ! Short function for Hulst-Curtis-Godson temperature factors for
3244 : ! computing effective H2O path
3245 : ! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
3246 : ! lines between 500 and 2820 cm^-1.
3247 : ! See cfa-www.harvard.edu/HITRAN
3248 : ! Isotopes of H2O: all
3249 : ! Line widths: air-broadened only (self set to 0)
3250 : ! Code for line strengths and widths: GENLN3
3251 : ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
3252 : ! Transmittance and Radiance Model, Version 3.0 Description
3253 : ! and Users Guide, NCAR/TN-367+STR, 147 pp.
3254 : !
3255 : ! Note: functions have been normalized by dividing by their values at
3256 : ! a path temperature of 160K
3257 : !
3258 : ! spectral intervals:
3259 : ! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
3260 : ! 2 = 800-1200 cm^-1
3261 : !
3262 : ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
3263 : ! 2nd edition, Oxford University Press, 1989.
3264 : ! Phi: function for H2O path
3265 : ! eq. 6.25, p. 228
3266 : !
3267 : real(r8),intent(in):: tpx ! path temperature
3268 : integer, intent(in):: iband ! band to process
3269 : real(r8) phi ! phi for given band
3270 : real(r8),parameter :: phi_r0(nbands) = (/ 9.60917711E-01_r8, -2.21031342E+01_r8/)
3271 : real(r8),parameter :: phi_r1(nbands) = (/ 4.86076751E-04_r8, 4.24062610E-01_r8/)
3272 : real(r8),parameter :: phi_r2(nbands) = (/-1.84806265E-06_r8, -2.95543415E-03_r8/)
3273 : real(r8),parameter :: phi_r3(nbands) = (/ 2.11239959E-09_r8, 7.52470896E-06_r8/)
3274 :
3275 2193998400 : phi = (((phi_r3(iband) * tpx) + phi_r2(iband)) * tpx + phi_r1(iband)) &
3276 2193998400 : * tpx + phi_r0(iband)
3277 2193998400 : end function phi
3278 :
3279 : !====================================================================================
3280 :
3281 2193998400 : function psi(tpx,iband)
3282 : !
3283 : ! History: First version for Hitran 1996 (C/H/E)
3284 : ! Current version for Hitran 2000 (C/LT/E)
3285 : ! Short function for Hulst-Curtis-Godson temperature factors for
3286 : ! computing effective H2O path
3287 : ! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
3288 : ! lines between 500 and 2820 cm^-1.
3289 : ! See cfa-www.harvard.edu/HITRAN
3290 : ! Isotopes of H2O: all
3291 : ! Line widths: air-broadened only (self set to 0)
3292 : ! Code for line strengths and widths: GENLN3
3293 : ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
3294 : ! Transmittance and Radiance Model, Version 3.0 Description
3295 : ! and Users Guide, NCAR/TN-367+STR, 147 pp.
3296 : !
3297 : ! Note: functions have been normalized by dividing by their values at
3298 : ! a path temperature of 160K
3299 : !
3300 : ! spectral intervals:
3301 : ! 1 = 0-800 cm^-1 and 1200-2200 cm^-1
3302 : ! 2 = 800-1200 cm^-1
3303 : !
3304 : ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
3305 : ! 2nd edition, Oxford University Press, 1989.
3306 : ! Psi: function for pressure along path
3307 : ! eq. 6.30, p. 228
3308 : !
3309 : real(r8),intent(in):: tpx ! path temperature
3310 : integer, intent(in):: iband ! band to process
3311 : real(r8) psi ! psi for given band
3312 : real(r8),parameter :: psi_r0(nbands) = (/ 5.65308452E-01_r8, -7.30087891E+01_r8/)
3313 : real(r8),parameter :: psi_r1(nbands) = (/ 4.07519005E-03_r8, 1.22199547E+00_r8/)
3314 : real(r8),parameter :: psi_r2(nbands) = (/-1.04347237E-05_r8, -7.12256227E-03_r8/)
3315 : real(r8),parameter :: psi_r3(nbands) = (/ 1.23765354E-08_r8, 1.47852825E-05_r8/)
3316 :
3317 2193998400 : psi = (((psi_r3(iband) * tpx) + psi_r2(iband)) * tpx + psi_r1(iband)) * tpx + psi_r0(iband)
3318 :
3319 2193998400 : end function psi
3320 :
3321 : !====================================================================================
3322 :
3323 47814624 : subroutine trcab(ncol , &
3324 : k1 ,k2 ,ucfc11 ,ucfc12 ,un2o0 , &
3325 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
3326 : uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
3327 : bch4 ,to3co2 ,pnm ,dw ,pnew , &
3328 : s2c ,uptype ,dplh2o ,abplnk1 ,tco2 , &
3329 : th2o ,to3 ,abstrc , &
3330 : aer_trn_ttl)
3331 : !-----------------------------------------------------------------------
3332 : !
3333 : ! Purpose:
3334 : ! Calculate absorptivity for non nearest layers for CH4, N2O, CFC11 and
3335 : ! CFC12.
3336 : !
3337 : ! Method:
3338 : ! See CCM3 description for equations.
3339 : !
3340 : ! Author: J. Kiehl
3341 : !
3342 : !-----------------------------------------------------------------------
3343 :
3344 : !------------------------------Arguments--------------------------------
3345 : !
3346 : ! Input arguments
3347 : !
3348 : integer, intent(in) :: ncol ! number of atmospheric columns
3349 : integer, intent(in) :: k1,k2 ! level indices
3350 : !
3351 : real(r8), intent(in) :: to3co2(pcols) ! pressure weighted temperature
3352 : real(r8), intent(in) :: pnm(pcols,pverp) ! interface pressures
3353 : real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
3354 : real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
3355 : real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
3356 : !
3357 : real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
3358 : real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
3359 : real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
3360 : real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
3361 : real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
3362 : !
3363 : real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
3364 : real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
3365 : real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
3366 : real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
3367 : real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
3368 : !
3369 : real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
3370 : real(r8), intent(in) :: dw(pcols) ! h2o path length
3371 : real(r8), intent(in) :: pnew(pcols) ! pressure
3372 : real(r8), intent(in) :: s2c(pcols,pverp) ! continuum path length
3373 : real(r8), intent(in) :: uptype(pcols,pverp) ! p-type h2o path length
3374 : !
3375 : real(r8), intent(in) :: dplh2o(pcols) ! p squared h2o path length
3376 : real(r8), intent(in) :: abplnk1(14,pcols,pverp) ! Planck factor
3377 : real(r8), intent(in) :: tco2(pcols) ! co2 transmission factor
3378 : real(r8), intent(in) :: th2o(pcols) ! h2o transmission factor
3379 : real(r8), intent(in) :: to3(pcols) ! o3 transmission factor
3380 :
3381 : real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands) ! aer trn.
3382 :
3383 : !
3384 : ! Output Arguments
3385 : !
3386 : real(r8), intent(out) :: abstrc(pcols) ! total trace gas absorptivity
3387 : !
3388 : !--------------------------Local Variables------------------------------
3389 : !
3390 : integer i,l ! loop counters
3391 :
3392 : real(r8) sqti(pcols) ! square root of mean temp
3393 : real(r8) du1 ! cfc11 path length
3394 : real(r8) du2 ! cfc12 path length
3395 : real(r8) acfc1 ! cfc11 absorptivity 798 cm-1
3396 : real(r8) acfc2 ! cfc11 absorptivity 846 cm-1
3397 : !
3398 : real(r8) acfc3 ! cfc11 absorptivity 933 cm-1
3399 : real(r8) acfc4 ! cfc11 absorptivity 1085 cm-1
3400 : real(r8) acfc5 ! cfc12 absorptivity 889 cm-1
3401 : real(r8) acfc6 ! cfc12 absorptivity 923 cm-1
3402 : real(r8) acfc7 ! cfc12 absorptivity 1102 cm-1
3403 : !
3404 : real(r8) acfc8 ! cfc12 absorptivity 1161 cm-1
3405 : real(r8) du01 ! n2o path length
3406 : real(r8) dbeta01 ! n2o pressure factor
3407 : real(r8) dbeta11 ! "
3408 : real(r8) an2o1 ! absorptivity of 1285 cm-1 n2o band
3409 : !
3410 : real(r8) du02 ! n2o path length
3411 : real(r8) dbeta02 ! n2o pressure factor
3412 : real(r8) an2o2 ! absorptivity of 589 cm-1 n2o band
3413 : real(r8) du03 ! n2o path length
3414 : real(r8) dbeta03 ! n2o pressure factor
3415 : !
3416 : real(r8) an2o3 ! absorptivity of 1168 cm-1 n2o band
3417 : real(r8) duch4 ! ch4 path length
3418 : real(r8) dbetac ! ch4 pressure factor
3419 : real(r8) ach4 ! absorptivity of 1306 cm-1 ch4 band
3420 : real(r8) du11 ! co2 path length
3421 : !
3422 : real(r8) du12 ! "
3423 : real(r8) du13 ! "
3424 : real(r8) dbetc1 ! co2 pressure factor
3425 : real(r8) dbetc2 ! co2 pressure factor
3426 : real(r8) aco21 ! absorptivity of 1064 cm-1 band
3427 : !
3428 : real(r8) du21 ! co2 path length
3429 : real(r8) du22 ! "
3430 : real(r8) du23 ! "
3431 : real(r8) aco22 ! absorptivity of 961 cm-1 band
3432 : real(r8) tt(pcols) ! temp. factor for h2o overlap factor
3433 : !
3434 : real(r8) psi1 ! "
3435 : real(r8) phi1 ! "
3436 : real(r8) p1 ! h2o overlap factor
3437 : real(r8) w1 ! "
3438 : real(r8) ds2c(pcols) ! continuum path length
3439 : !
3440 : real(r8) duptyp(pcols) ! p-type path length
3441 : real(r8) tw(pcols,6) ! h2o transmission factor
3442 : ! real(r8) g1(6) ! "
3443 : ! real(r8) g2(6) ! "
3444 : ! real(r8) g3(6) ! "
3445 : !
3446 : ! real(r8) g4(6) ! "
3447 : ! real(r8) ab(6) ! h2o temp. factor
3448 : ! real(r8) bb(6) ! "
3449 : ! real(r8) abp(6) ! "
3450 : ! real(r8) bbp(6) ! "
3451 : !
3452 : real(r8) tcfc3 ! transmission for cfc11 band
3453 : real(r8) tcfc4 ! transmission for cfc11 band
3454 : real(r8) tcfc6 ! transmission for cfc12 band
3455 : real(r8) tcfc7 ! transmission for cfc12 band
3456 : real(r8) tcfc8 ! transmission for cfc12 band
3457 : !
3458 : real(r8) tlw ! h2o transmission
3459 : real(r8) tch4 ! ch4 transmission
3460 : !
3461 : !--------------------------Data Statements------------------------------
3462 : !
3463 : ! data g1 /0.0468556_r8,0.0397454_r8,0.0407664_r8,0.0304380_r8,0.0540398_r8,0.0321962_r8/
3464 : ! data g2 /14.4832_r8,4.30242_r8,5.23523_r8,3.25342_r8,0.698935_r8,16.5599_r8/
3465 : ! data g3 /26.1898_r8,18.4476_r8,15.3633_r8,12.1927_r8,9.14992_r8,8.07092_r8/
3466 : ! data g4 /0.0261782_r8,0.0369516_r8,0.0307266_r8,0.0243854_r8,0.0182932_r8,0.0161418_r8/
3467 : ! data ab /3.0857e-2_r8,2.3524e-2_r8,1.7310e-2_r8,2.6661e-2_r8,2.8074e-2_r8,2.2915e-2_r8/
3468 : ! data bb /-1.3512e-4_r8,-6.8320e-5_r8,-3.2609e-5_r8,-1.0228e-5_r8,-9.5743e-5_r8,-1.0304e-4_r8/
3469 : ! data abp/2.9129e-2_r8,2.4101e-2_r8,1.9821e-2_r8,2.6904e-2_r8,2.9458e-2_r8,1.9892e-2_r8/
3470 : ! data bbp/-1.3139e-4_r8,-5.5688e-5_r8,-4.6380e-5_r8,-8.0362e-5_r8,-1.0115e-4_r8,-8.8061e-5_r8/
3471 : !
3472 : !--------------------------Statement Functions--------------------------
3473 : !
3474 : real(r8) func, u, b
3475 : func(u,b) = u/sqrt(4.0_r8 + u*(1.0_r8 + 1.0_r8 / b))
3476 : !
3477 : !------------------------------------------------------------------------
3478 : !
3479 798393024 : do i = 1,ncol
3480 750578400 : sqti(i) = sqrt(to3co2(i))
3481 : !
3482 : ! h2o transmission
3483 : !
3484 750578400 : tt(i) = abs(to3co2(i) - 250.0_r8)
3485 750578400 : ds2c(i) = abs(s2c(i,k1) - s2c(i,k2))
3486 798393024 : duptyp(i) = abs(uptype(i,k1) - uptype(i,k2))
3487 : end do
3488 : !
3489 334702368 : do l = 1,6
3490 4838172768 : do i = 1,ncol
3491 4503470400 : psi1 = exp(abp(l)*tt(i) + bbp(l)*tt(i)*tt(i))
3492 4503470400 : phi1 = exp(ab(l)*tt(i) + bb(l)*tt(i)*tt(i))
3493 4503470400 : p1 = pnew(i)*(psi1/phi1)/sslp
3494 4503470400 : w1 = dw(i)*phi1
3495 : tw(i,l) = exp(-g1(l)*p1*(sqrt(1.0_r8 + g2(l)*(w1/p1)) - 1.0_r8) - &
3496 4790358144 : g3(l)*ds2c(i)-g4(l)*duptyp(i))
3497 : end do
3498 : end do
3499 : !
3500 798393024 : do i=1,ncol
3501 2251735200 : tw(i,1)=tw(i,1)*(0.7_r8*aer_trn_ttl(i,k1,k2,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
3502 2251735200 : 0.3_r8*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000))
3503 750578400 : tw(i,2)=tw(i,2)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
3504 750578400 : tw(i,3)=tw(i,3)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
3505 750578400 : tw(i,4)=tw(i,4)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
3506 750578400 : tw(i,5)=tw(i,5)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
3507 798393024 : tw(i,6)=tw(i,6)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
3508 : end do ! end loop over lon
3509 798393024 : do i = 1,ncol
3510 750578400 : du1 = abs(ucfc11(i,k1) - ucfc11(i,k2))
3511 750578400 : du2 = abs(ucfc12(i,k1) - ucfc12(i,k2))
3512 : !
3513 : ! cfc transmissions
3514 : !
3515 750578400 : tcfc3 = exp(-175.005_r8*du1)
3516 750578400 : tcfc4 = exp(-1202.18_r8*du1)
3517 750578400 : tcfc6 = exp(-5786.73_r8*du2)
3518 750578400 : tcfc7 = exp(-2873.51_r8*du2)
3519 750578400 : tcfc8 = exp(-2085.59_r8*du2)
3520 : !
3521 : ! Absorptivity for CFC11 bands
3522 : !
3523 750578400 : acfc1 = 50.0_r8*(1.0_r8 - exp(-54.09_r8*du1))*tw(i,1)*abplnk1(7,i,k2)
3524 750578400 : acfc2 = 60.0_r8*(1.0_r8 - exp(-5130.03_r8*du1))*tw(i,2)*abplnk1(8,i,k2)
3525 750578400 : acfc3 = 60.0_r8*(1.0_r8 - tcfc3)*tw(i,4)*tcfc6*abplnk1(9,i,k2)
3526 750578400 : acfc4 = 100.0_r8*(1.0_r8 - tcfc4)*tw(i,5)*abplnk1(10,i,k2)
3527 : !
3528 : ! Absorptivity for CFC12 bands
3529 : !
3530 750578400 : acfc5 = 45.0_r8*(1.0_r8 - exp(-1272.35_r8*du2))*tw(i,3)*abplnk1(11,i,k2)
3531 750578400 : acfc6 = 50.0_r8*(1.0_r8 - tcfc6)* tw(i,4) * abplnk1(12,i,k2)
3532 750578400 : acfc7 = 80.0_r8*(1.0_r8 - tcfc7)* tw(i,5) * tcfc4*abplnk1(13,i,k2)
3533 750578400 : acfc8 = 70.0_r8*(1.0_r8 - tcfc8)* tw(i,6) * abplnk1(14,i,k2)
3534 : !
3535 : ! Emissivity for CH4 band 1306 cm-1
3536 : !
3537 750578400 : tlw = exp(-1.0_r8*sqrt(dplh2o(i)))
3538 750578400 : tlw=tlw*aer_trn_ttl(i,k1,k2,idx_LW_1200_2000)
3539 750578400 : duch4 = abs(uch4(i,k1) - uch4(i,k2))
3540 750578400 : dbetac = abs(bch4(i,k1) - bch4(i,k2))/duch4
3541 750578400 : ach4 = 6.00444_r8*sqti(i)*log(1.0_r8 + func(duch4,dbetac))*tlw*abplnk1(3,i,k2)
3542 750578400 : tch4 = 1.0_r8/(1.0_r8 + 0.02_r8*func(duch4,dbetac))
3543 : !
3544 : ! Absorptivity for N2O bands
3545 : !
3546 750578400 : du01 = abs(un2o0(i,k1) - un2o0(i,k2))
3547 750578400 : du11 = abs(un2o1(i,k1) - un2o1(i,k2))
3548 750578400 : dbeta01 = abs(bn2o0(i,k1) - bn2o0(i,k2))/du01
3549 750578400 : dbeta11 = abs(bn2o1(i,k1) - bn2o1(i,k2))/du11
3550 : !
3551 : ! 1285 cm-1 band
3552 : !
3553 : an2o1 = 2.35558_r8*sqti(i)*log(1.0_r8 + func(du01,dbeta01) &
3554 750578400 : + func(du11,dbeta11))*tlw*tch4*abplnk1(4,i,k2)
3555 750578400 : du02 = 0.100090_r8*du01
3556 750578400 : du12 = 0.0992746_r8*du11
3557 750578400 : dbeta02 = 0.964282_r8*dbeta01
3558 : !
3559 : ! 589 cm-1 band
3560 : !
3561 : an2o2 = 2.65581_r8*sqti(i)*log(1.0_r8 + func(du02,dbeta02) + &
3562 750578400 : func(du12,dbeta02))*th2o(i)*tco2(i)*abplnk1(5,i,k2)
3563 750578400 : du03 = 0.0333767_r8*du01
3564 750578400 : dbeta03 = 0.982143_r8*dbeta01
3565 : !
3566 : ! 1168 cm-1 band
3567 : !
3568 : an2o3 = 2.54034_r8*sqti(i)*log(1.0_r8 + func(du03,dbeta03))* &
3569 750578400 : tw(i,6)*tcfc8*abplnk1(6,i,k2)
3570 : !
3571 : ! Emissivity for 1064 cm-1 band of CO2
3572 : !
3573 750578400 : du11 = abs(uco211(i,k1) - uco211(i,k2))
3574 750578400 : du12 = abs(uco212(i,k1) - uco212(i,k2))
3575 750578400 : du13 = abs(uco213(i,k1) - uco213(i,k2))
3576 750578400 : dbetc1 = 2.97558_r8*abs(pnm(i,k1) + pnm(i,k2))/(2.0_r8*sslp*sqti(i))
3577 750578400 : dbetc2 = 2.0_r8*dbetc1
3578 : aco21 = 3.7571_r8*sqti(i)*log(1.0_r8 + func(du11,dbetc1) &
3579 : + func(du12,dbetc2) + func(du13,dbetc2)) &
3580 750578400 : *to3(i)*tw(i,5)*tcfc4*tcfc7*abplnk1(2,i,k2)
3581 : !
3582 : ! Emissivity for 961 cm-1 band
3583 : !
3584 750578400 : du21 = abs(uco221(i,k1) - uco221(i,k2))
3585 750578400 : du22 = abs(uco222(i,k1) - uco222(i,k2))
3586 750578400 : du23 = abs(uco223(i,k1) - uco223(i,k2))
3587 : aco22 = 3.8443_r8*sqti(i)*log(1.0_r8 + func(du21,dbetc1) &
3588 : + func(du22,dbetc1) + func(du23,dbetc2)) &
3589 750578400 : *tw(i,4)*tcfc3*tcfc6*abplnk1(1,i,k2)
3590 : !
3591 : ! total trace gas absorptivity
3592 : !
3593 : abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
3594 : acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
3595 798393024 : aco21 + aco22
3596 : end do
3597 :
3598 47814624 : end subroutine trcab
3599 :
3600 : !====================================================================================
3601 :
3602 7083648 : subroutine trcabn(ncol , &
3603 : k2 ,kn ,ucfc11 ,ucfc12 ,un2o0 , &
3604 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
3605 : uco221 ,uco222 ,uco223 ,tbar ,bplnk , &
3606 : winpl ,pinpl ,tco2 ,th2o ,to3 , &
3607 : uptype ,dw ,s2c ,up2 ,pnew , &
3608 : abstrc ,uinpl , &
3609 : aer_trn_ngh)
3610 : !-----------------------------------------------------------------------
3611 : !
3612 : ! Purpose:
3613 : ! Calculate nearest layer absorptivity due to CH4, N2O, CFC11 and CFC12
3614 : !
3615 : ! Method:
3616 : ! Equations in CCM3 description
3617 : !
3618 : ! Author: J. Kiehl
3619 : !
3620 : !-----------------------------------------------------------------------
3621 :
3622 : !------------------------------Arguments--------------------------------
3623 : !
3624 : ! Input arguments
3625 : !
3626 : integer, intent(in) :: ncol ! number of atmospheric columns
3627 : integer, intent(in) :: k2 ! level index
3628 : integer, intent(in) :: kn ! level index
3629 : !
3630 : real(r8), intent(in) :: tbar(pcols,4) ! pressure weighted temperature
3631 : real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
3632 : real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
3633 : real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
3634 : real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
3635 : !
3636 : real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
3637 : real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
3638 : real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
3639 : real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
3640 : real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
3641 : !
3642 : real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
3643 : real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
3644 : real(r8), intent(in) :: bplnk(14,pcols,4) ! weighted Planck fnc. for absorptivity
3645 : real(r8), intent(in) :: winpl(pcols,4) ! fractional path length
3646 : real(r8), intent(in) :: pinpl(pcols,4) ! pressure factor for subdivided layer
3647 : !
3648 : real(r8), intent(in) :: tco2(pcols) ! co2 transmission
3649 : real(r8), intent(in) :: th2o(pcols) ! h2o transmission
3650 : real(r8), intent(in) :: to3(pcols) ! o3 transmission
3651 : real(r8), intent(in) :: dw(pcols) ! h2o path length
3652 : real(r8), intent(in) :: pnew(pcols) ! pressure factor
3653 : !
3654 : real(r8), intent(in) :: s2c(pcols,pverp) ! h2o continuum factor
3655 : real(r8), intent(in) :: uptype(pcols,pverp) ! p-type path length
3656 : real(r8), intent(in) :: up2(pcols) ! p squared path length
3657 : real(r8), intent(in) :: uinpl(pcols,4) ! Nearest layer subdivision factor
3658 : real(r8), intent(in) :: aer_trn_ngh(pcols,nlwbands)
3659 : ! [fraction] Total transmission between
3660 : ! nearest neighbor sub-levels
3661 : !
3662 : ! Output Arguments
3663 : !
3664 : real(r8), intent(out) :: abstrc(pcols) ! total trace gas absorptivity
3665 :
3666 : !
3667 : !--------------------------Local Variables------------------------------
3668 : !
3669 : integer i,l ! loop counters
3670 : !
3671 : real(r8) sqti(pcols) ! square root of mean temp
3672 : real(r8) rsqti(pcols) ! reciprocal of sqti
3673 : real(r8) du1 ! cfc11 path length
3674 : real(r8) du2 ! cfc12 path length
3675 : real(r8) acfc1 ! absorptivity of cfc11 798 cm-1 band
3676 : !
3677 : real(r8) acfc2 ! absorptivity of cfc11 846 cm-1 band
3678 : real(r8) acfc3 ! absorptivity of cfc11 933 cm-1 band
3679 : real(r8) acfc4 ! absorptivity of cfc11 1085 cm-1 band
3680 : real(r8) acfc5 ! absorptivity of cfc11 889 cm-1 band
3681 : real(r8) acfc6 ! absorptivity of cfc11 923 cm-1 band
3682 : !
3683 : real(r8) acfc7 ! absorptivity of cfc11 1102 cm-1 band
3684 : real(r8) acfc8 ! absorptivity of cfc11 1161 cm-1 band
3685 : real(r8) du01 ! n2o path length
3686 : real(r8) dbeta01 ! n2o pressure factors
3687 : real(r8) dbeta11 ! "
3688 : !
3689 : real(r8) an2o1 ! absorptivity of the 1285 cm-1 n2o band
3690 : real(r8) du02 ! n2o path length
3691 : real(r8) dbeta02 ! n2o pressure factor
3692 : real(r8) an2o2 ! absorptivity of the 589 cm-1 n2o band
3693 : real(r8) du03 ! n2o path length
3694 : !
3695 : real(r8) dbeta03 ! n2o pressure factor
3696 : real(r8) an2o3 ! absorptivity of the 1168 cm-1 n2o band
3697 : real(r8) duch4 ! ch4 path length
3698 : real(r8) dbetac ! ch4 pressure factor
3699 : real(r8) ach4 ! absorptivity of the 1306 cm-1 ch4 band
3700 : !
3701 : real(r8) du11 ! co2 path length
3702 : real(r8) du12 ! "
3703 : real(r8) du13 ! "
3704 : real(r8) dbetc1 ! co2 pressure factor
3705 : real(r8) dbetc2 ! co2 pressure factor
3706 : !
3707 : real(r8) aco21 ! absorptivity of the 1064 cm-1 co2 band
3708 : real(r8) du21 ! co2 path length
3709 : real(r8) du22 ! "
3710 : real(r8) du23 ! "
3711 : real(r8) aco22 ! absorptivity of the 961 cm-1 co2 band
3712 : !
3713 : real(r8) tt(pcols) ! temp. factor for h2o overlap
3714 : real(r8) psi1 ! "
3715 : real(r8) phi1 ! "
3716 : real(r8) p1 ! factor for h2o overlap
3717 : real(r8) w1 ! "
3718 : !
3719 : real(r8) ds2c(pcols) ! continuum path length
3720 : real(r8) duptyp(pcols) ! p-type path length
3721 : real(r8) tw(pcols,6) ! h2o transmission overlap
3722 : ! real(r8) g1(6) ! h2o overlap factor
3723 : ! real(r8) g2(6) ! "
3724 : !
3725 : ! real(r8) g3(6) ! "
3726 : ! real(r8) g4(6) ! "
3727 : ! real(r8) ab(6) ! h2o temp. factor
3728 : ! real(r8) bb(6) ! "
3729 : ! real(r8) abp(6) ! "
3730 : !
3731 : ! real(r8) bbp(6) ! "
3732 : real(r8) tcfc3 ! transmission of cfc11 band
3733 : real(r8) tcfc4 ! transmission of cfc11 band
3734 : real(r8) tcfc6 ! transmission of cfc12 band
3735 : real(r8) tcfc7 ! "
3736 : !
3737 : real(r8) tcfc8 ! "
3738 : real(r8) tlw ! h2o transmission
3739 : real(r8) tch4 ! ch4 transmission
3740 : !
3741 : !--------------------------Data Statements------------------------------
3742 : !
3743 : ! data g1 /0.0468556_r8,0.0397454_r8,0.0407664_r8,0.0304380_r8,0.0540398_r8,0.0321962_r8/
3744 : ! data g2 /14.4832_r8,4.30242_r8,5.23523_r8,3.25342_r8,0.698935_r8,16.5599_r8/
3745 : ! data g3 /26.1898_r8,18.4476_r8,15.3633_r8,12.1927_r8,9.14992_r8,8.07092_r8/
3746 : ! data g4 /0.0261782_r8,0.0369516_r8,0.0307266_r8,0.0243854_r8,0.0182932_r8,0.0161418_r8/
3747 : ! data ab /3.0857e-2_r8,2.3524e-2_r8,1.7310e-2_r8,2.6661e-2_r8,2.8074e-2_r8,2.2915e-2_r8/
3748 : ! data bb /-1.3512e-4_r8,-6.8320e-5_r8,-3.2609e-5_r8,-1.0228e-5_r8,-9.5743e-5_r8,-1.0304e-4_r8/
3749 : ! data abp/2.9129e-2_r8,2.4101e-2_r8,1.9821e-2_r8,2.6904e-2_r8,2.9458e-2_r8,1.9892e-2_r8/
3750 : ! data bbp/-1.3139e-4_r8,-5.5688e-5_r8,-4.6380e-5_r8,-8.0362e-5_r8,-1.0115e-4_r8,-8.8061e-5_r8/
3751 : !
3752 : !--------------------------Statement Functions--------------------------
3753 : !
3754 : real(r8) func, u, b
3755 : func(u,b) = u/sqrt(4.0_r8 + u*(1.0_r8 + 1.0_r8 / b))
3756 : !
3757 : !------------------------------------------------------------------
3758 : !
3759 118280448 : do i = 1,ncol
3760 111196800 : sqti(i) = sqrt(tbar(i,kn))
3761 111196800 : rsqti(i) = 1._r8 / sqti(i)
3762 : !
3763 : ! h2o transmission
3764 : !
3765 111196800 : tt(i) = abs(tbar(i,kn) - 250.0_r8)
3766 111196800 : ds2c(i) = abs(s2c(i,k2+1) - s2c(i,k2))*uinpl(i,kn)
3767 118280448 : duptyp(i) = abs(uptype(i,k2+1) - uptype(i,k2))*uinpl(i,kn)
3768 : end do
3769 : !
3770 49585536 : do l = 1,6
3771 716766336 : do i = 1,ncol
3772 667180800 : psi1 = exp(abp(l)*tt(i)+bbp(l)*tt(i)*tt(i))
3773 667180800 : phi1 = exp(ab(l)*tt(i)+bb(l)*tt(i)*tt(i))
3774 667180800 : p1 = pnew(i) * (psi1/phi1) / sslp
3775 667180800 : w1 = dw(i) * winpl(i,kn) * phi1
3776 : tw(i,l) = exp(- g1(l)*p1*(sqrt(1.0_r8+g2(l)*(w1/p1))-1.0_r8) &
3777 709682688 : - g3(l)*ds2c(i)-g4(l)*duptyp(i))
3778 : end do
3779 : end do
3780 : !
3781 118280448 : do i=1,ncol
3782 111196800 : tw(i,1)=tw(i,1)*(0.7_r8*aer_trn_ngh(i,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
3783 111196800 : 0.3_r8*aer_trn_ngh(i,idx_LW_0800_1000))
3784 111196800 : tw(i,2)=tw(i,2)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
3785 111196800 : tw(i,3)=tw(i,3)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
3786 111196800 : tw(i,4)=tw(i,4)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
3787 111196800 : tw(i,5)=tw(i,5)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
3788 118280448 : tw(i,6)=tw(i,6)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
3789 : end do ! end loop over lon
3790 :
3791 118280448 : do i = 1,ncol
3792 : !
3793 111196800 : du1 = abs(ucfc11(i,k2+1) - ucfc11(i,k2)) * winpl(i,kn)
3794 111196800 : du2 = abs(ucfc12(i,k2+1) - ucfc12(i,k2)) * winpl(i,kn)
3795 : !
3796 : ! cfc transmissions
3797 : !
3798 111196800 : tcfc3 = exp(-175.005_r8*du1)
3799 111196800 : tcfc4 = exp(-1202.18_r8*du1)
3800 111196800 : tcfc6 = exp(-5786.73_r8*du2)
3801 111196800 : tcfc7 = exp(-2873.51_r8*du2)
3802 111196800 : tcfc8 = exp(-2085.59_r8*du2)
3803 : !
3804 : ! Absorptivity for CFC11 bands
3805 : !
3806 111196800 : acfc1 = 50.0_r8*(1.0_r8 - exp(-54.09_r8*du1)) * tw(i,1)*bplnk(7,i,kn)
3807 111196800 : acfc2 = 60.0_r8*(1.0_r8 - exp(-5130.03_r8*du1))*tw(i,2)*bplnk(8,i,kn)
3808 111196800 : acfc3 = 60.0_r8*(1.0_r8 - tcfc3)*tw(i,4)*tcfc6 * bplnk(9,i,kn)
3809 111196800 : acfc4 = 100.0_r8*(1.0_r8 - tcfc4)* tw(i,5) * bplnk(10,i,kn)
3810 : !
3811 : ! Absorptivity for CFC12 bands
3812 : !
3813 111196800 : acfc5 = 45.0_r8*(1.0_r8 - exp(-1272.35_r8*du2))*tw(i,3)*bplnk(11,i,kn)
3814 111196800 : acfc6 = 50.0_r8*(1.0_r8 - tcfc6)*tw(i,4)*bplnk(12,i,kn)
3815 111196800 : acfc7 = 80.0_r8*(1.0_r8 - tcfc7)* tw(i,5)*tcfc4 *bplnk(13,i,kn)
3816 111196800 : acfc8 = 70.0_r8*(1.0_r8 - tcfc8)*tw(i,6)*bplnk(14,i,kn)
3817 : !
3818 : ! Absorptivity for CH4 band 1306 cm-1
3819 : !
3820 111196800 : tlw = exp(-1.0_r8*sqrt(up2(i)))
3821 111196800 : tlw=tlw*aer_trn_ngh(i,idx_LW_1200_2000)
3822 111196800 : duch4 = abs(uch4(i,k2+1) - uch4(i,k2)) * winpl(i,kn)
3823 111196800 : dbetac = 2.94449_r8 * pinpl(i,kn) * rsqti(i) / sslp
3824 111196800 : ach4 = 6.00444_r8*sqti(i)*log(1.0_r8 + func(duch4,dbetac)) * tlw * bplnk(3,i,kn)
3825 111196800 : tch4 = 1.0_r8/(1.0_r8 + 0.02_r8*func(duch4,dbetac))
3826 : !
3827 : ! Absorptivity for N2O bands
3828 : !
3829 111196800 : du01 = abs(un2o0(i,k2+1) - un2o0(i,k2)) * winpl(i,kn)
3830 111196800 : du11 = abs(un2o1(i,k2+1) - un2o1(i,k2)) * winpl(i,kn)
3831 111196800 : dbeta01 = 19.399_r8 * pinpl(i,kn) * rsqti(i) / sslp
3832 111196800 : dbeta11 = dbeta01
3833 : !
3834 : ! 1285 cm-1 band
3835 : !
3836 : an2o1 = 2.35558_r8*sqti(i)*log(1.0_r8 + func(du01,dbeta01) &
3837 111196800 : + func(du11,dbeta11)) * tlw * tch4 * bplnk(4,i,kn)
3838 111196800 : du02 = 0.100090_r8*du01
3839 111196800 : du12 = 0.0992746_r8*du11
3840 111196800 : dbeta02 = 0.964282_r8*dbeta01
3841 : !
3842 : ! 589 cm-1 band
3843 : !
3844 : an2o2 = 2.65581_r8*sqti(i)*log(1.0_r8 + func(du02,dbeta02) &
3845 111196800 : + func(du12,dbeta02)) * tco2(i) * th2o(i) * bplnk(5,i,kn)
3846 111196800 : du03 = 0.0333767_r8*du01
3847 111196800 : dbeta03 = 0.982143_r8*dbeta01
3848 : !
3849 : ! 1168 cm-1 band
3850 : !
3851 : an2o3 = 2.54034_r8*sqti(i)*log(1.0_r8 + func(du03,dbeta03)) * &
3852 111196800 : tw(i,6) * tcfc8 * bplnk(6,i,kn)
3853 : !
3854 : ! Absorptivity for 1064 cm-1 band of CO2
3855 : !
3856 111196800 : du11 = abs(uco211(i,k2+1) - uco211(i,k2)) * winpl(i,kn)
3857 111196800 : du12 = abs(uco212(i,k2+1) - uco212(i,k2)) * winpl(i,kn)
3858 111196800 : du13 = abs(uco213(i,k2+1) - uco213(i,k2)) * winpl(i,kn)
3859 111196800 : dbetc1 = 2.97558_r8 * pinpl(i,kn) * rsqti(i) / sslp
3860 111196800 : dbetc2 = 2.0_r8 * dbetc1
3861 : aco21 = 3.7571_r8*sqti(i)*log(1.0_r8 + func(du11,dbetc1) &
3862 : + func(du12,dbetc2) + func(du13,dbetc2)) &
3863 111196800 : * to3(i) * tw(i,5) * tcfc4 * tcfc7 * bplnk(2,i,kn)
3864 : !
3865 : ! Absorptivity for 961 cm-1 band of co2
3866 : !
3867 111196800 : du21 = abs(uco221(i,k2+1) - uco221(i,k2)) * winpl(i,kn)
3868 111196800 : du22 = abs(uco222(i,k2+1) - uco222(i,k2)) * winpl(i,kn)
3869 111196800 : du23 = abs(uco223(i,k2+1) - uco223(i,k2)) * winpl(i,kn)
3870 : aco22 = 3.8443_r8*sqti(i)*log(1.0_r8 + func(du21,dbetc1) &
3871 : + func(du22,dbetc1) + func(du23,dbetc2)) &
3872 111196800 : * tw(i,4) * tcfc3 * tcfc6 * bplnk(1,i,kn)
3873 : !
3874 : ! total trace gas absorptivity
3875 : !
3876 : abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
3877 : acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
3878 118280448 : aco21 + aco22
3879 : end do
3880 :
3881 7083648 : end subroutine trcabn
3882 :
3883 : !====================================================================================
3884 :
3885 1839024 : subroutine trcems(ncol , &
3886 : k ,co2t ,pnm ,ucfc11 ,ucfc12 , &
3887 : un2o0 ,un2o1 ,bn2o0 ,bn2o1 ,uch4 , &
3888 : bch4 ,uco211 ,uco212 ,uco213 ,uco221 , &
3889 : uco222 ,uco223 ,uptype ,w ,s2c , &
3890 : up2 ,emplnk ,th2o ,tco2 ,to3 , &
3891 : emstrc , &
3892 : aer_trn_ttl)
3893 : !-----------------------------------------------------------------------
3894 : !
3895 : ! Purpose:
3896 : ! Calculate emissivity for CH4, N2O, CFC11 and CFC12 bands.
3897 : !
3898 : ! Method:
3899 : ! See CCM3 Description for equations.
3900 : !
3901 : ! Author: J. Kiehl
3902 : !
3903 : !-----------------------------------------------------------------------
3904 :
3905 : !------------------------------Arguments--------------------------------
3906 : !
3907 : ! Input arguments
3908 : !
3909 : integer, intent(in) :: ncol ! number of atmospheric columns
3910 :
3911 : real(r8), intent(in) :: co2t(pcols,pverp) ! pressure weighted temperature
3912 : real(r8), intent(in) :: pnm(pcols,pverp) ! interface pressure
3913 : real(r8), intent(in) :: ucfc11(pcols,pverp) ! CFC11 path length
3914 : real(r8), intent(in) :: ucfc12(pcols,pverp) ! CFC12 path length
3915 : real(r8), intent(in) :: un2o0(pcols,pverp) ! N2O path length
3916 : !
3917 : real(r8), intent(in) :: un2o1(pcols,pverp) ! N2O path length (hot band)
3918 : real(r8), intent(in) :: uch4(pcols,pverp) ! CH4 path length
3919 : real(r8), intent(in) :: uco211(pcols,pverp) ! CO2 9.4 micron band path length
3920 : real(r8), intent(in) :: uco212(pcols,pverp) ! CO2 9.4 micron band path length
3921 : real(r8), intent(in) :: uco213(pcols,pverp) ! CO2 9.4 micron band path length
3922 : !
3923 : real(r8), intent(in) :: uco221(pcols,pverp) ! CO2 10.4 micron band path length
3924 : real(r8), intent(in) :: uco222(pcols,pverp) ! CO2 10.4 micron band path length
3925 : real(r8), intent(in) :: uco223(pcols,pverp) ! CO2 10.4 micron band path length
3926 : real(r8), intent(in) :: uptype(pcols,pverp) ! continuum path length
3927 : real(r8), intent(in) :: bn2o0(pcols,pverp) ! pressure factor for n2o
3928 : !
3929 : real(r8), intent(in) :: bn2o1(pcols,pverp) ! pressure factor for n2o
3930 : real(r8), intent(in) :: bch4(pcols,pverp) ! pressure factor for ch4
3931 : real(r8), intent(in) :: emplnk(14,pcols) ! emissivity Planck factor
3932 : real(r8), intent(in) :: th2o(pcols) ! water vapor overlap factor
3933 : real(r8), intent(in) :: tco2(pcols) ! co2 overlap factor
3934 : !
3935 : real(r8), intent(in) :: to3(pcols) ! o3 overlap factor
3936 : real(r8), intent(in) :: s2c(pcols,pverp) ! h2o continuum path length
3937 : real(r8), intent(in) :: w(pcols,pverp) ! h2o path length
3938 : real(r8), intent(in) :: up2(pcols) ! pressure squared h2o path length
3939 : !
3940 : integer, intent(in) :: k ! level index
3941 :
3942 : real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands) ! aer trn.
3943 :
3944 : !
3945 : ! Output Arguments
3946 : !
3947 : real(r8), intent(out) :: emstrc(pcols,pverp) ! total trace gas emissivity
3948 :
3949 : !
3950 : !--------------------------Local Variables------------------------------
3951 : !
3952 : integer i,l ! loop counters
3953 : !
3954 : real(r8) sqti(pcols) ! square root of mean temp
3955 : real(r8) ecfc1 ! emissivity of cfc11 798 cm-1 band
3956 : real(r8) ecfc2 ! " " " 846 cm-1 band
3957 : real(r8) ecfc3 ! " " " 933 cm-1 band
3958 : real(r8) ecfc4 ! " " " 1085 cm-1 band
3959 : !
3960 : real(r8) ecfc5 ! " " cfc12 889 cm-1 band
3961 : real(r8) ecfc6 ! " " " 923 cm-1 band
3962 : real(r8) ecfc7 ! " " " 1102 cm-1 band
3963 : real(r8) ecfc8 ! " " " 1161 cm-1 band
3964 : real(r8) u01 ! n2o path length
3965 : !
3966 : real(r8) u11 ! n2o path length
3967 : real(r8) beta01 ! n2o pressure factor
3968 : real(r8) beta11 ! n2o pressure factor
3969 : real(r8) en2o1 ! emissivity of the 1285 cm-1 N2O band
3970 : real(r8) u02 ! n2o path length
3971 : !
3972 : real(r8) u12 ! n2o path length
3973 : real(r8) beta02 ! n2o pressure factor
3974 : real(r8) en2o2 ! emissivity of the 589 cm-1 N2O band
3975 : real(r8) u03 ! n2o path length
3976 : real(r8) beta03 ! n2o pressure factor
3977 : !
3978 : real(r8) en2o3 ! emissivity of the 1168 cm-1 N2O band
3979 : real(r8) betac ! ch4 pressure factor
3980 : real(r8) ech4 ! emissivity of 1306 cm-1 CH4 band
3981 : real(r8) betac1 ! co2 pressure factor
3982 : real(r8) betac2 ! co2 pressure factor
3983 : !
3984 : real(r8) eco21 ! emissivity of 1064 cm-1 CO2 band
3985 : real(r8) eco22 ! emissivity of 961 cm-1 CO2 band
3986 : real(r8) tt(pcols) ! temp. factor for h2o overlap factor
3987 : real(r8) psi1 ! narrow band h2o temp. factor
3988 : real(r8) phi1 ! "
3989 : !
3990 : real(r8) p1 ! h2o line overlap factor
3991 : real(r8) w1 ! "
3992 : real(r8) tw(pcols,6) ! h2o transmission overlap
3993 : ! real(r8) g1(6) ! h2o overlap factor
3994 : ! real(r8) g2(6) ! "
3995 : !
3996 : ! real(r8) g3(6) ! "
3997 : ! real(r8) g4(6) ! "
3998 : ! real(r8) ab(6) ! "
3999 : ! real(r8) bb(6) ! "
4000 : ! real(r8) abp(6) ! "
4001 : !
4002 : ! real(r8) bbp(6) ! "
4003 : real(r8) tcfc3 ! transmission for cfc11 band
4004 : real(r8) tcfc4 ! "
4005 : real(r8) tcfc6 ! transmission for cfc12 band
4006 : real(r8) tcfc7 ! "
4007 : !
4008 : real(r8) tcfc8 ! "
4009 : real(r8) tlw ! h2o overlap factor
4010 : real(r8) tch4 ! ch4 overlap factor
4011 : !
4012 : !--------------------------Data Statements------------------------------
4013 : !
4014 : ! data g1 /0.0468556_r8,0.0397454_r8,0.0407664_r8,0.0304380_r8,0.0540398_r8,0.0321962_r8/
4015 : ! data g2 /14.4832_r8,4.30242_r8,5.23523_r8,3.25342_r8,0.698935_r8,16.5599_r8/
4016 : ! data g3 /26.1898_r8,18.4476_r8,15.3633_r8,12.1927_r8,9.14992_r8,8.07092_r8/
4017 : ! data g4 /0.0261782_r8,0.0369516_r8,0.0307266_r8,0.0243854_r8,0.0182932_r8,0.0161418_r8/
4018 : ! data ab /3.0857e-2_r8,2.3524e-2_r8,1.7310e-2_r8,2.6661e-2_r8,2.8074e-2_r8,2.2915e-2_r8/
4019 : ! data bb /-1.3512e-4_r8,-6.8320e-5_r8,-3.2609e-5_r8,-1.0228e-5_r8,-9.5743e-5_r8,-1.0304e-4_r8/
4020 : ! data abp/2.9129e-2_r8,2.4101e-2_r8,1.9821e-2_r8,2.6904e-2_r8,2.9458e-2_r8,1.9892e-2_r8/
4021 : ! data bbp/-1.3139e-4_r8,-5.5688e-5_r8,-4.6380e-5_r8,-8.0362e-5_r8,-1.0115e-4_r8,-8.8061e-5_r8/
4022 : !
4023 : !--------------------------Statement Functions--------------------------
4024 : !
4025 : real(r8) func, u, b
4026 : func(u,b) = u/sqrt(4.0_r8 + u*(1.0_r8 + 1.0_r8 / b))
4027 : !
4028 : !-----------------------------------------------------------------------
4029 : !
4030 30707424 : do i = 1,ncol
4031 28868400 : sqti(i) = sqrt(co2t(i,k))
4032 : !
4033 : ! Transmission for h2o
4034 : !
4035 30707424 : tt(i) = abs(co2t(i,k) - 250.0_r8)
4036 : end do
4037 : !
4038 12873168 : do l = 1,6
4039 186083568 : do i = 1,ncol
4040 173210400 : psi1 = exp(abp(l)*tt(i)+bbp(l)*tt(i)*tt(i))
4041 173210400 : phi1 = exp(ab(l)*tt(i)+bb(l)*tt(i)*tt(i))
4042 173210400 : p1 = pnm(i,k) * (psi1/phi1) / sslp
4043 173210400 : w1 = w(i,k) * phi1
4044 : tw(i,l) = exp(- g1(l)*p1*(sqrt(1.0_r8+g2(l)*(w1/p1))-1.0_r8) &
4045 184244544 : - g3(l)*s2c(i,k)-g4(l)*uptype(i,k))
4046 : end do
4047 : end do
4048 :
4049 : ! Overlap H2O tranmission with STRAER continuum in 6 trace gas
4050 : ! subbands
4051 :
4052 30707424 : do i=1,ncol
4053 57736800 : tw(i,1)=tw(i,1)*(0.7_r8*aer_trn_ttl(i,k,1,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
4054 57736800 : 0.3_r8*aer_trn_ttl(i,k,1,idx_LW_0800_1000))
4055 28868400 : tw(i,2)=tw(i,2)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
4056 28868400 : tw(i,3)=tw(i,3)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
4057 28868400 : tw(i,4)=tw(i,4)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
4058 28868400 : tw(i,5)=tw(i,5)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
4059 30707424 : tw(i,6)=tw(i,6)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
4060 : end do ! end loop over lon
4061 : !
4062 30707424 : do i = 1,ncol
4063 : !
4064 : ! transmission due to cfc bands
4065 : !
4066 28868400 : tcfc3 = exp(-175.005_r8*ucfc11(i,k))
4067 28868400 : tcfc4 = exp(-1202.18_r8*ucfc11(i,k))
4068 28868400 : tcfc6 = exp(-5786.73_r8*ucfc12(i,k))
4069 28868400 : tcfc7 = exp(-2873.51_r8*ucfc12(i,k))
4070 28868400 : tcfc8 = exp(-2085.59_r8*ucfc12(i,k))
4071 : !
4072 : ! Emissivity for CFC11 bands
4073 : !
4074 28868400 : ecfc1 = 50.0_r8*(1.0_r8 - exp(-54.09_r8*ucfc11(i,k))) * tw(i,1) * emplnk(7,i)
4075 28868400 : ecfc2 = 60.0_r8*(1.0_r8 - exp(-5130.03_r8*ucfc11(i,k)))* tw(i,2) * emplnk(8,i)
4076 28868400 : ecfc3 = 60.0_r8*(1.0_r8 - tcfc3)*tw(i,4)*tcfc6*emplnk(9,i)
4077 28868400 : ecfc4 = 100.0_r8*(1.0_r8 - tcfc4)*tw(i,5)*emplnk(10,i)
4078 : !
4079 : ! Emissivity for CFC12 bands
4080 : !
4081 28868400 : ecfc5 = 45.0_r8*(1.0_r8 - exp(-1272.35_r8*ucfc12(i,k)))*tw(i,3)*emplnk(11,i)
4082 28868400 : ecfc6 = 50.0_r8*(1.0_r8 - tcfc6)*tw(i,4)*emplnk(12,i)
4083 28868400 : ecfc7 = 80.0_r8*(1.0_r8 - tcfc7)*tw(i,5)* tcfc4 * emplnk(13,i)
4084 28868400 : ecfc8 = 70.0_r8*(1.0_r8 - tcfc8)*tw(i,6) * emplnk(14,i)
4085 : !
4086 : ! Emissivity for CH4 band 1306 cm-1
4087 : !
4088 28868400 : tlw = exp(-1.0_r8*sqrt(up2(i)))
4089 :
4090 : ! Overlap H2O vibration rotation band with STRAER continuum
4091 : ! for CH4 1306 cm-1 and N2O 1285 cm-1 bands
4092 :
4093 28868400 : tlw=tlw*aer_trn_ttl(i,k,1,idx_LW_1200_2000)
4094 28868400 : betac = bch4(i,k)/uch4(i,k)
4095 28868400 : ech4 = 6.00444_r8*sqti(i)*log(1.0_r8 + func(uch4(i,k),betac)) *tlw * emplnk(3,i)
4096 28868400 : tch4 = 1.0_r8/(1.0_r8 + 0.02_r8*func(uch4(i,k),betac))
4097 : !
4098 : ! Emissivity for N2O bands
4099 : !
4100 28868400 : u01 = un2o0(i,k)
4101 28868400 : u11 = un2o1(i,k)
4102 28868400 : beta01 = bn2o0(i,k)/un2o0(i,k)
4103 28868400 : beta11 = bn2o1(i,k)/un2o1(i,k)
4104 : !
4105 : ! 1285 cm-1 band
4106 : !
4107 : en2o1 = 2.35558_r8*sqti(i)*log(1.0_r8 + func(u01,beta01) + &
4108 28868400 : func(u11,beta11))*tlw*tch4*emplnk(4,i)
4109 28868400 : u02 = 0.100090_r8*u01
4110 28868400 : u12 = 0.0992746_r8*u11
4111 28868400 : beta02 = 0.964282_r8*beta01
4112 : !
4113 : ! 589 cm-1 band
4114 : !
4115 : en2o2 = 2.65581_r8*sqti(i)*log(1.0_r8 + func(u02,beta02) + &
4116 28868400 : func(u12,beta02)) * tco2(i) * th2o(i) * emplnk(5,i)
4117 28868400 : u03 = 0.0333767_r8*u01
4118 28868400 : beta03 = 0.982143_r8*beta01
4119 : !
4120 : ! 1168 cm-1 band
4121 : !
4122 : en2o3 = 2.54034_r8*sqti(i)*log(1.0_r8 + func(u03,beta03)) * &
4123 28868400 : tw(i,6) * tcfc8 * emplnk(6,i)
4124 : !
4125 : ! Emissivity for 1064 cm-1 band of CO2
4126 : !
4127 28868400 : betac1 = 2.97558_r8*pnm(i,k) / (sslp*sqti(i))
4128 28868400 : betac2 = 2.0_r8 * betac1
4129 : eco21 = 3.7571_r8*sqti(i)*log(1.0_r8 + func(uco211(i,k),betac1) &
4130 : + func(uco212(i,k),betac2) + func(uco213(i,k),betac2)) &
4131 28868400 : * to3(i) * tw(i,5) * tcfc4 * tcfc7 * emplnk(2,i)
4132 : !
4133 : ! Emissivity for 961 cm-1 band
4134 : !
4135 : eco22 = 3.8443_r8*sqti(i)*log(1.0_r8 + func(uco221(i,k),betac1) &
4136 : + func(uco222(i,k),betac1) + func(uco223(i,k),betac2)) &
4137 28868400 : * tw(i,4) * tcfc3 * tcfc6 * emplnk(1,i)
4138 : !
4139 : ! total trace gas emissivity
4140 : !
4141 : emstrc(i,k) = ecfc1 + ecfc2 + ecfc3 + ecfc4 + ecfc5 +ecfc6 + &
4142 : ecfc7 + ecfc8 + en2o1 + en2o2 + en2o3 + ech4 + &
4143 30707424 : eco21 + eco22
4144 : end do
4145 :
4146 1839024 : end subroutine trcems
4147 :
4148 : !====================================================================================
4149 :
4150 68112 : subroutine trcplk(ncol , &
4151 : tint ,tlayr ,tplnke ,emplnk ,abplnk1 , &
4152 : abplnk2 )
4153 : !-----------------------------------------------------------------------
4154 : !
4155 : ! Purpose:
4156 : ! Calculate Planck factors for absorptivity and emissivity of
4157 : ! CH4, N2O, CFC11 and CFC12
4158 : !
4159 : ! Method:
4160 : ! Planck function and derivative evaluated at the band center.
4161 : !
4162 : ! Author: J. Kiehl
4163 : !
4164 : !------------------------------Arguments--------------------------------
4165 : !
4166 : ! Input arguments
4167 : !
4168 : integer, intent(in) :: ncol ! number of atmospheric columns
4169 :
4170 : real(r8), intent(in) :: tint(pcols,pverp) ! interface temperatures
4171 : real(r8), intent(in) :: tlayr(pcols,pverp) ! k-1 level temperatures
4172 : real(r8), intent(in) :: tplnke(pcols) ! Top Layer temperature
4173 : !
4174 : ! output arguments
4175 : !
4176 : real(r8), intent(out) :: emplnk(14,pcols) ! emissivity Planck factor
4177 : real(r8), intent(out) :: abplnk1(14,pcols,pverp) ! non-nearest layer Plack factor
4178 : real(r8), intent(out) :: abplnk2(14,pcols,pverp) ! nearest layer factor
4179 :
4180 : !
4181 : !--------------------------Local Variables------------------------------
4182 : !
4183 : integer wvl ! wavelength index
4184 : integer i,k ! loop counters
4185 : !
4186 : real(r8) f1(14) ! Planck function factor
4187 : real(r8) f2(14) ! "
4188 : real(r8) f3(14) ! "
4189 : !
4190 : !--------------------------Data Statements------------------------------
4191 : !
4192 : data f1 /5.85713e8_r8,7.94950e8_r8,1.47009e9_r8,1.40031e9_r8,1.34853e8_r8, &
4193 : 1.05158e9_r8,3.35370e8_r8,3.99601e8_r8,5.35994e8_r8,8.42955e8_r8, &
4194 : 4.63682e8_r8,5.18944e8_r8,8.83202e8_r8,1.03279e9_r8/
4195 : data f2 /2.02493e11_r8,3.04286e11_r8,6.90698e11_r8,6.47333e11_r8, &
4196 : 2.85744e10_r8,4.41862e11_r8,9.62780e10_r8,1.21618e11_r8, &
4197 : 1.79905e11_r8,3.29029e11_r8,1.48294e11_r8,1.72315e11_r8, &
4198 : 3.50140e11_r8,4.31364e11_r8/
4199 : data f3 /1383.0_r8,1531.0_r8,1879.0_r8,1849.0_r8,848.0_r8,1681.0_r8, &
4200 : 1148.0_r8,1217.0_r8,1343.0_r8,1561.0_r8,1279.0_r8,1328.0_r8, &
4201 : 1586.0_r8,1671.0_r8/
4202 : !
4203 : !-----------------------------------------------------------------------
4204 : !
4205 : ! Calculate emissivity Planck factor
4206 : !
4207 1021680 : do wvl = 1,14
4208 15990480 : do i = 1,ncol
4209 15922368 : emplnk(wvl,i) = f1(wvl)/(tplnke(i)**4.0_r8*(exp(f3(wvl)/tplnke(i))-1.0_r8))
4210 : end do
4211 : end do
4212 : !
4213 : ! Calculate absorptivity Planck factor for tint and tlayr temperatures
4214 : !
4215 1021680 : do wvl = 1,14
4216 26768016 : do k = ntoplw, pverp
4217 430857504 : do i = 1, ncol
4218 : !
4219 : ! non-nearlest layer function
4220 : !
4221 1212472800 : abplnk1(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tint(i,k))) &
4222 1212472800 : /(tint(i,k)**5.0_r8*(exp(f3(wvl)/tint(i,k))-1.0_r8)**2.0_r8)
4223 : !
4224 : ! nearest layer function
4225 : !
4226 : abplnk2(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tlayr(i,k))) &
4227 429903936 : /(tlayr(i,k)**5.0_r8*(exp(f3(wvl)/tlayr(i,k))-1.0_r8)**2.0_r8)
4228 : end do
4229 : end do
4230 : end do
4231 :
4232 68112 : end subroutine trcplk
4233 :
4234 :
4235 : !====================================================================================
4236 :
4237 : end module radae
|
