Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rrtmg_lw_rad.nomcica.f90 @ 4318

Last change on this file since 4318 was 1585, checked in by maronga, 10 years ago
Added support for RRTMG radiation code
File size: 48.8 KB

Line
1	! path: $Source: /storm/rc1/cvsroot/rc/rrtmg_lw/src/rrtmg_lw_rad.nomcica.f90,v $
2	! author: $Author: miacono $
3	! revision: $Revision: 1.12 $
4	! created: $Date: 2011/04/08 20:25:01 $
5	!
6
7	module rrtmg_lw_rad
8
9	! --------------------------------------------------------------------------
10	! \| \|
11	! \| Copyright 2002-2009, Atmospheric & Environmental Research, Inc. (AER). \|
12	! \| This software may be used, copied, or redistributed as long as it is \|
13	! \| not sold and this copyright notice is reproduced on each copy made. \|
14	! \| This model is provided as is without any express or implied warranties. \|
15	! \| (http://www.rtweb.aer.com/) \|
16	! \| \|
17	! --------------------------------------------------------------------------
18	!
19	! ****************************************************************************
20	! * *
21	! * RRTMG_LW *
22	! * *
23	! * *
24	! * *
25	! * a rapid radiative transfer model *
26	! * for the longwave region *
27	! * for application to general circulation models *
28	! * *
29	! * *
30	! * Atmospheric and Environmental Research, Inc. *
31	! * 131 Hartwell Avenue *
32	! * Lexington, MA 02421 *
33	! * *
34	! * *
35	! * Eli J. Mlawer *
36	! * Jennifer S. Delamere *
37	! * Michael J. Iacono *
38	! * Shepard A. Clough *
39	! * *
40	! * *
41	! * *
42	! * *
43	! * *
44	! * *
45	! * email: emlawer@aer.com *
46	! * email: jdelamer@aer.com *
47	! * email: miacono@aer.com *
48	! * *
49	! * The authors wish to acknowledge the contributions of the *
50	! * following people: Steven J. Taubman, Karen Cady-Pereira, *
51	! * Patrick D. Brown, Ronald E. Farren, Luke Chen, Robert Bergstrom. *
52	! * *
53	! ****************************************************************************
54
55	! -------- Modules --------
56	use parkind, only : im => kind_im, rb => kind_rb
57	use rrlw_vsn
58	use rrtmg_lw_cldprop, only: cldprop
59	! *** Move the required call to rrtmg_lw_ini below and the following
60	! use association to the GCM initialization area ***
61	! use rrtmg_lw_init, only: rrtmg_lw_ini
62	use rrtmg_lw_rtrn, only: rtrn
63	use rrtmg_lw_rtrnmr, only: rtrnmr
64	use rrtmg_lw_setcoef, only: setcoef
65	use rrtmg_lw_taumol, only: taumol
66
67	implicit none
68
69	! public interfaces/functions/subroutines
70	public :: rrtmg_lw, inatm
71
72	!------------------------------------------------------------------
73	contains
74	!------------------------------------------------------------------
75
76	!------------------------------------------------------------------
77	! Public subroutines
78	!------------------------------------------------------------------
79
80	subroutine rrtmg_lw &
81	(ncol ,nlay ,icld ,idrv , &
82	play ,plev ,tlay ,tlev ,tsfc , &
83	h2ovmr ,o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr, &
84	cfc11vmr,cfc12vmr,cfc22vmr,ccl4vmr ,emis , &
85	inflglw ,iceflglw,liqflglw,cldfr , &
86	taucld ,cicewp ,cliqwp ,reice ,reliq , &
87	tauaer , &
88	uflx ,dflx ,hr ,uflxc ,dflxc, hrc, &
89	duflx_dt,duflxc_dt )
90
91	! -------- Description --------
92
93	! This program is the driver subroutine for RRTMG_LW, the AER LW radiation
94	! model for application to GCMs, that has been adapted from RRTM_LW for
95	! improved efficiency.
96	!
97	! NOTE: The call to RRTMG_LW_INI should be moved to the GCM initialization
98	! area, since this has to be called only once.
99	!
100	! This routine:
101	! a) calls INATM to read in the atmospheric profile from GCM;
102	! all layering in RRTMG is ordered from surface to toa.
103	! b) calls CLDPROP to set cloud optical depth based on input
104	! cloud properties
105	! c) calls SETCOEF to calculate various quantities needed for
106	! the radiative transfer algorithm
107	! d) calls TAUMOL to calculate gaseous optical depths for each
108	! of the 16 spectral bands
109	! e) calls RTRNMR (for both clear and cloudy profiles) to perform the
110	! radiative transfer calculation with a maximum-random cloud
111	! overlap method, or calls RTRN to use random cloud overlap.
112	! f) passes the necessary fluxes and cooling rates back to GCM
113	!
114	! Two modes of operation are possible:
115	! The mode is chosen by using either rrtmg_lw.nomcica.f90 (to not use
116	! McICA) or rrtmg_lw.f90 (to use McICA) to interface with a GCM.
117	!
118	! 1) Standard, single forward model calculation (imca = 0)
119	! 2) Monte Carlo Independent Column Approximation (McICA, Pincus et al.,
120	! JC, 2003) method is applied to the forward model calculation (imca = 1)
121	!
122	! Two methods of cloud property input are possible:
123	! Cloud properties can be input in one of two ways (controlled by input
124	! flags inflglw, iceflglw, and liqflglw; see text file rrtmg_lw_instructions
125	! and subroutine rrtmg_lw_cldprop.f90 for further details):
126	!
127	! 1) Input cloud fraction and cloud optical depth directly (inflglw = 0)
128	! 2) Input cloud fraction and cloud physical properties (inflglw = 1 or 2);
129	! cloud optical properties are calculated by cldprop or cldprmc based
130	! on input settings of iceflglw and liqflglw. Ice particle size provided
131	! must be appropriately defined for the ice parameterization selected.
132	!
133	! One method of aerosol property input is possible:
134	! Aerosol properties can be input in only one way (controlled by input
135	! flag iaer; see text file rrtmg_lw_instructions for further details):
136	!
137	! 1) Input aerosol optical depth directly by layer and spectral band (iaer=10);
138	! band average optical depth at the mid-point of each spectral band.
139	! RRTMG_LW currently treats only aerosol absorption;
140	! scattering capability is not presently available.
141	!
142	! The optional calculation of the change in upward flux as a function of surface
143	! temperature is available (controlled by input flag idrv). This can be utilized
144	! to approximate adjustments to the upward flux profile caused only by a change in
145	! surface temperature between full radiation calls. This feature uses the pre-
146	! calculated derivative of the Planck function with respect to surface temperature.
147	!
148	! 1) Normal forward calculation for the input profile (idrv=0)
149	! 2) Normal forward calculation with optional calculation of the change
150	! in upward flux as a function of surface temperature for clear sky
151	! and total sky flux. Flux partial derivatives are provided in arrays
152	! duflx_dt and duflxc_dt for total and clear sky. (idrv=1)
153	!
154	!
155	! ------- Modifications -------
156	!
157	! This version of RRTMG_LW has been modified from RRTM_LW to use a reduced
158	! set of g-points for application to GCMs.
159	!
160	!-- Original version (derived from RRTM_LW), reduction of g-points, other
161	! revisions for use with GCMs.
162	! 1999: M. J. Iacono, AER, Inc.
163	!-- Adapted for use with NCAR/CAM.
164	! May 2004: M. J. Iacono, AER, Inc.
165	!-- Conversion to F90 formatting for consistency with rrtmg_sw.
166	! Feb 2007: M. J. Iacono, AER, Inc.
167	!-- Modifications to formatting to use assumed-shape arrays.
168	! Aug 2007: M. J. Iacono, AER, Inc.
169	!-- Modified to add longwave aerosol absorption.
170	! Apr 2008: M. J. Iacono, AER, Inc.
171	!-- Added capability to calculate derivative of upward flux wrt surface temperature.
172	! Nov 2009: M. J. Iacono, E. J. Mlawer, AER, Inc.
173
174	! --------- Modules ----------
175
176	use parrrtm, only : nbndlw, ngptlw, maxxsec, mxmol
177	use rrlw_con, only: fluxfac, heatfac, oneminus, pi
178	use rrlw_wvn, only: ng, ngb, nspa, nspb, wavenum1, wavenum2, delwave
179
180	! ------- Declarations -------
181
182	! ----- Input -----
183	! Note: All volume mixing ratios are in dimensionless units of mole fraction obtained
184	! by scaling mass mixing ratio (g/g) with the appropriate molecular weights (g/mol)
185	integer(kind=im), intent(in) :: ncol ! Number of horizontal columns
186	integer(kind=im), intent(in) :: nlay ! Number of model layers
187	integer(kind=im), intent(inout) :: icld ! Cloud overlap method
188	! 0: Clear only
189	! 1: Random
190	! 2: Maximum/random
191	! 3: Maximum
192	integer(kind=im), intent(in) :: idrv ! Flag for calculation of dFdT, the change
193	! in upward flux as a function of
194	! surface temperature [0=off, 1=on]
195	! 0: Normal forward calculation
196	! 1: Normal forward calculation with
197	! duflx_dt and duflxc_dt output
198
199	real(kind=rb), intent(in) :: play(:,:) ! Layer pressures (hPa, mb)
200	! Dimensions: (ncol,nlay)
201	real(kind=rb), intent(in) :: plev(:,:) ! Interface pressures (hPa, mb)
202	! Dimensions: (ncol,nlay+1)
203	real(kind=rb), intent(in) :: tlay(:,:) ! Layer temperatures (K)
204	! Dimensions: (ncol,nlay)
205	real(kind=rb), intent(in) :: tlev(:,:) ! Interface temperatures (K)
206	! Dimensions: (ncol,nlay+1)
207	real(kind=rb), intent(in) :: tsfc(:) ! Surface temperature (K)
208	! Dimensions: (ncol)
209	real(kind=rb), intent(in) :: h2ovmr(:,:) ! H2O volume mixing ratio
210	! Dimensions: (ncol,nlay)
211	real(kind=rb), intent(in) :: o3vmr(:,:) ! O3 volume mixing ratio
212	! Dimensions: (ncol,nlay)
213	real(kind=rb), intent(in) :: co2vmr(:,:) ! CO2 volume mixing ratio
214	! Dimensions: (ncol,nlay)
215	real(kind=rb), intent(in) :: ch4vmr(:,:) ! Methane volume mixing ratio
216	! Dimensions: (ncol,nlay)
217	real(kind=rb), intent(in) :: n2ovmr(:,:) ! Nitrous oxide volume mixing ratio
218	! Dimensions: (ncol,nlay)
219	real(kind=rb), intent(in) :: o2vmr(:,:) ! Oxygen volume mixing ratio
220	! Dimensions: (ncol,nlay)
221	real(kind=rb), intent(in) :: cfc11vmr(:,:) ! CFC11 volume mixing ratio
222	! Dimensions: (ncol,nlay)
223	real(kind=rb), intent(in) :: cfc12vmr(:,:) ! CFC12 volume mixing ratio
224	! Dimensions: (ncol,nlay)
225	real(kind=rb), intent(in) :: cfc22vmr(:,:) ! CFC22 volume mixing ratio
226	! Dimensions: (ncol,nlay)
227	real(kind=rb), intent(in) :: ccl4vmr(:,:) ! CCL4 volume mixing ratio
228	! Dimensions: (ncol,nlay)
229	real(kind=rb), intent(in) :: emis(:,:) ! Surface emissivity
230	! Dimensions: (ncol,nbndlw)
231
232	integer(kind=im), intent(in) :: inflglw ! Flag for cloud optical properties
233	integer(kind=im), intent(in) :: iceflglw ! Flag for ice particle specification
234	integer(kind=im), intent(in) :: liqflglw ! Flag for liquid droplet specification
235
236	real(kind=rb), intent(in) :: cldfr(:,:) ! Cloud fraction
237	! Dimensions: (ncol,nlay)
238	real(kind=rb), intent(in) :: cicewp(:,:) ! Cloud ice water path (g/m2)
239	! Dimensions: (ncol,nlay)
240	real(kind=rb), intent(in) :: cliqwp(:,:) ! Cloud liquid water path (g/m2)
241	! Dimensions: (ncol,nlay)
242	real(kind=rb), intent(in) :: reice(:,:) ! Cloud ice particle effective size (microns)
243	! Dimensions: (ncol,nlay)
244	! specific definition of reice depends on setting of iceflglw:
245	! iceflglw = 0: ice effective radius, r_ec, (Ebert and Curry, 1992),
246	! r_ec must be >= 10.0 microns
247	! iceflglw = 1: ice effective radius, r_ec, (Ebert and Curry, 1992),
248	! r_ec range is limited to 13.0 to 130.0 microns
249	! iceflglw = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996)
250	! r_k range is limited to 5.0 to 131.0 microns
251	! iceflglw = 3: generalized effective size, dge, (Fu, 1996),
252	! dge range is limited to 5.0 to 140.0 microns
253	! [dge = 1.0315 * r_ec]
254	real(kind=rb), intent(in) :: reliq(:,:) ! Cloud water drop effective radius (microns)
255	! Dimensions: (ncol,nlay)
256	real(kind=rb), intent(in) :: taucld(:,:,:) ! In-cloud optical depth
257	! Dimensions: (nbndlw,ncol,nlay)
258	! real(kind=rb), intent(in) :: ssacld(:,:,:) ! In-cloud single scattering albedo
259	! Dimensions: (nbndlw,ncol,nlay)
260	! for future expansion
261	! (lw scattering not yet available)
262	! real(kind=rb), intent(in) :: asmcld(:,:,:) ! In-cloud asymmetry parameter
263	! Dimensions: (nbndlw,ncol,nlay)
264	! for future expansion
265	! (lw scattering not yet available)
266	real(kind=rb), intent(in) :: tauaer(:,:,:) ! aerosol optical depth
267	! Dimensions: (ncol,nlay,nbndlw)
268	! real(kind=rb), intent(in) :: ssaaer(:,:,:) ! aerosol single scattering albedo
269	! Dimensions: (ncol,nlay,nbndlw)
270	! for future expansion
271	! (lw aerosols/scattering not yet available)
272	! real(kind=rb), intent(in) :: asmaer(:,:,:) ! aerosol asymmetry parameter
273	! Dimensions: (ncol,nlay,nbndlw)
274	! for future expansion
275	! (lw aerosols/scattering not yet available)
276
277
278	! ----- Output -----
279
280	real(kind=rb), intent(out) :: uflx(:,:) ! Total sky longwave upward flux (W/m2)
281	! Dimensions: (ncol,nlay+1)
282	real(kind=rb), intent(out) :: dflx(:,:) ! Total sky longwave downward flux (W/m2)
283	! Dimensions: (ncol,nlay+1)
284	real(kind=rb), intent(out) :: hr(:,:) ! Total sky longwave radiative heating rate (K/d)
285	! Dimensions: (ncol,nlay)
286	real(kind=rb), intent(out) :: uflxc(:,:) ! Clear sky longwave upward flux (W/m2)
287	! Dimensions: (ncol,nlay+1)
288	real(kind=rb), intent(out) :: dflxc(:,:) ! Clear sky longwave downward flux (W/m2)
289	! Dimensions: (ncol,nlay+1)
290	real(kind=rb), intent(out) :: hrc(:,:) ! Clear sky longwave radiative heating rate (K/d)
291	! Dimensions: (ncol,nlay)
292
293	! ----- Optional Output -----
294	real(kind=rb), intent(out), optional :: duflx_dt(:,:)
295	! change in upward longwave flux (w/m2/k)
296	! with respect to surface temperature
297	! Dimensions: (ncol,nlay+1)
298	real(kind=rb), intent(out), optional :: duflxc_dt(:,:)
299	! change in clear sky upward longwave flux (w/m2/k)
300	! with respect to surface temperature
301	! Dimensions: (ncol,nlay+1)
302
303	! ----- Local -----
304
305	! Control
306	integer(kind=im) :: nlayers ! total number of layers
307	integer(kind=im) :: istart ! beginning band of calculation
308	integer(kind=im) :: iend ! ending band of calculation
309	integer(kind=im) :: iout ! output option flag (inactive)
310	integer(kind=im) :: iaer ! aerosol option flag
311	integer(kind=im) :: iplon ! column loop index
312	integer(kind=im) :: imca ! flag for mcica [0=off, 1=on]
313	integer(kind=im) :: k ! layer loop index
314	integer(kind=im) :: ig ! g-point loop index
315
316	! Atmosphere
317	real(kind=rb) :: pavel(nlay+1) ! layer pressures (mb)
318	real(kind=rb) :: tavel(nlay+1) ! layer temperatures (K)
319	real(kind=rb) :: pz(0:nlay+1) ! level (interface) pressures (hPa, mb)
320	real(kind=rb) :: tz(0:nlay+1) ! level (interface) temperatures (K)
321	real(kind=rb) :: tbound ! surface temperature (K)
322	real(kind=rb) :: coldry(nlay+1) ! dry air column density (mol/cm2)
323	real(kind=rb) :: wbrodl(nlay+1) ! broadening gas column density (mol/cm2)
324	real(kind=rb) :: wkl(mxmol,nlay+1) ! molecular amounts (mol/cm-2)
325	real(kind=rb) :: wx(maxxsec,nlay+1) ! cross-section amounts (mol/cm-2)
326	real(kind=rb) :: pwvcm ! precipitable water vapor (cm)
327	real(kind=rb) :: semiss(nbndlw) ! lw surface emissivity
328	real(kind=rb) :: fracs(nlay+1,ngptlw) !
329	real(kind=rb) :: taug(nlay+1,ngptlw) ! gaseous optical depths
330	real(kind=rb) :: taut(nlay+1,ngptlw) ! gaseous + aerosol optical depths
331
332	real(kind=rb) :: taua(nlay+1,nbndlw) ! aerosol optical depth
333	! real(kind=rb) :: ssaa(nlay+1,nbndlw) ! aerosol single scattering albedo
334	! for future expansion
335	! (lw aerosols/scattering not yet available)
336	! real(kind=rb) :: asma(nlay+1,nbndlw) ! aerosol asymmetry parameter
337	! for future expansion
338	! (lw aerosols/scattering not yet available)
339
340	! Atmosphere - setcoef
341	integer(kind=im) :: laytrop ! tropopause layer index
342	integer(kind=im) :: jp(nlay+1) ! lookup table index
343	integer(kind=im) :: jt(nlay+1) ! lookup table index
344	integer(kind=im) :: jt1(nlay+1) ! lookup table index
345	real(kind=rb) :: planklay(nlay+1,nbndlw)!
346	real(kind=rb) :: planklev(0:nlay+1,nbndlw)!
347	real(kind=rb) :: plankbnd(nbndlw) !
348	real(kind=rb) :: dplankbnd_dt(nbndlw) !
349
350	real(kind=rb) :: colh2o(nlay+1) ! column amount (h2o)
351	real(kind=rb) :: colco2(nlay+1) ! column amount (co2)
352	real(kind=rb) :: colo3(nlay+1) ! column amount (o3)
353	real(kind=rb) :: coln2o(nlay+1) ! column amount (n2o)
354	real(kind=rb) :: colco(nlay+1) ! column amount (co)
355	real(kind=rb) :: colch4(nlay+1) ! column amount (ch4)
356	real(kind=rb) :: colo2(nlay+1) ! column amount (o2)
357	real(kind=rb) :: colbrd(nlay+1) ! column amount (broadening gases)
358
359	integer(kind=im) :: indself(nlay+1)
360	integer(kind=im) :: indfor(nlay+1)
361	real(kind=rb) :: selffac(nlay+1)
362	real(kind=rb) :: selffrac(nlay+1)
363	real(kind=rb) :: forfac(nlay+1)
364	real(kind=rb) :: forfrac(nlay+1)
365
366	integer(kind=im) :: indminor(nlay+1)
367	real(kind=rb) :: minorfrac(nlay+1)
368	real(kind=rb) :: scaleminor(nlay+1)
369	real(kind=rb) :: scaleminorn2(nlay+1)
370
371	real(kind=rb) :: & !
372	fac00(nlay+1), fac01(nlay+1), &
373	fac10(nlay+1), fac11(nlay+1)
374	real(kind=rb) :: & !
375	rat_h2oco2(nlay+1),rat_h2oco2_1(nlay+1), &
376	rat_h2oo3(nlay+1),rat_h2oo3_1(nlay+1), &
377	rat_h2on2o(nlay+1),rat_h2on2o_1(nlay+1), &
378	rat_h2och4(nlay+1),rat_h2och4_1(nlay+1), &
379	rat_n2oco2(nlay+1),rat_n2oco2_1(nlay+1), &
380	rat_o3co2(nlay+1),rat_o3co2_1(nlay+1)
381
382	! Atmosphere/clouds - cldprop
383	integer(kind=im) :: ncbands ! number of cloud spectral bands
384	integer(kind=im) :: inflag ! flag for cloud property method
385	integer(kind=im) :: iceflag ! flag for ice cloud properties
386	integer(kind=im) :: liqflag ! flag for liquid cloud properties
387
388	real(kind=rb) :: cldfrac(nlay+1) ! layer cloud fraction
389	real(kind=rb) :: tauc(nbndlw,nlay+1) ! in-cloud optical depth
390	! real(kind=rb) :: ssac(nbndlw,nlay+1) ! in-cloud single scattering albedo
391	! for future expansion
392	! (lw scattering not yet available)
393	! real(kind=rb) :: asmc(nbndlw,nlay+1) ! in-cloud asymmetry parameter
394	! for future expansion
395	! (lw scattering not yet available)
396	real(kind=rb) :: ciwp(nlay+1) ! cloud ice water path
397	real(kind=rb) :: clwp(nlay+1) ! cloud liquid water path
398	real(kind=rb) :: rel(nlay+1) ! cloud liquid particle effective radius (microns)
399	real(kind=rb) :: rei(nlay+1) ! cloud ice particle effective size (microns)
400	real(kind=rb) :: taucloud(nlay+1,nbndlw)! layer in-cloud optical depth
401
402	! Output
403	real(kind=rb) :: totuflux(0:nlay+1) ! upward longwave flux (w/m2)
404	real(kind=rb) :: totdflux(0:nlay+1) ! downward longwave flux (w/m2)
405	real(kind=rb) :: fnet(0:nlay+1) ! net longwave flux (w/m2)
406	real(kind=rb) :: htr(0:nlay+1) ! longwave heating rate (k/day)
407	real(kind=rb) :: totuclfl(0:nlay+1) ! clear sky upward longwave flux (w/m2)
408	real(kind=rb) :: totdclfl(0:nlay+1) ! clear sky downward longwave flux (w/m2)
409	real(kind=rb) :: fnetc(0:nlay+1) ! clear sky net longwave flux (w/m2)
410	real(kind=rb) :: htrc(0:nlay+1) ! clear sky longwave heating rate (k/day)
411	real(kind=rb) :: dtotuflux_dt(0:nlay+1) ! change in upward longwave flux (w/m2/k)
412	! with respect to surface temperature
413	real(kind=rb) :: dtotuclfl_dt(0:nlay+1) ! change in clear sky upward longwave flux (w/m2/k)
414	! with respect to surface temperature
415
416	!
417	! Initializations
418
419	oneminus = 1._rb - 1.e-6_rb
420	pi = 2._rb*asin(1._rb)
421	fluxfac = pi * 2.e4_rb ! orig: fluxfac = pi * 2.d4
422	istart = 1
423	iend = 16
424	iout = 0
425
426	! Set imca to select calculation type:
427	! imca = 0, use standard forward model calculation
428	! imca = 1, use McICA for Monte Carlo treatment of sub-grid cloud variability
429
430	! * This version does not use McICA (imca = 0) *
431
432	! Set default icld to select of clear or cloud calculation and cloud overlap method
433	! icld = 0, clear only
434	! icld = 1, with clouds using random cloud overlap
435	! icld = 2, with clouds using maximum/random cloud overlap
436	! icld = 3, with clouds using maximum cloud overlap (McICA only)
437	if (icld.lt.0.or.icld.gt.3) icld = 2
438
439	! Set iaer to select aerosol option
440	! iaer = 0, no aerosols
441	! icld = 10, input total aerosol optical depth (tauaer) directly
442	iaer = 10
443
444	! Call model and data initialization, compute lookup tables, perform
445	! reduction of g-points from 256 to 140 for input absorption coefficient
446	! data and other arrays.
447	!
448	! In a GCM this call should be placed in the model initialization
449	! area, since this has to be called only once.
450	! call rrtmg_lw_ini(cpdair)
451
452	! This is the main longitude/column loop within RRTMG.
453	do iplon = 1, ncol
454
455	! Prepare atmospheric profile from GCM for use in RRTMG, and define
456	! other input parameters.
457
458	call inatm (iplon, nlay, icld, iaer, &
459	play, plev, tlay, tlev, tsfc, h2ovmr, &
460	o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr, cfc11vmr, cfc12vmr, &
461	cfc22vmr, ccl4vmr, emis, inflglw, iceflglw, liqflglw, &
462	cldfr, taucld, cicewp, cliqwp, reice, reliq, tauaer, &
463	nlayers, pavel, pz, tavel, tz, tbound, semiss, coldry, &
464	wkl, wbrodl, wx, pwvcm, inflag, iceflag, liqflag, &
465	cldfrac, tauc, ciwp, clwp, rei, rel, taua)
466
467	! For cloudy atmosphere, use cldprop to set cloud optical properties based on
468	! input cloud physical properties. Select method based on choices described
469	! in cldprop. Cloud fraction, water path, liquid droplet and ice particle
470	! effective radius must be passed into cldprop. Cloud fraction and cloud
471	! optical depth are transferred to rrtmg_lw arrays in cldprop.
472
473	call cldprop(nlayers, inflag, iceflag, liqflag, cldfrac, tauc, &
474	ciwp, clwp, rei, rel, ncbands, taucloud)
475
476	! Calculate information needed by the radiative transfer routine
477	! that is specific to this atmosphere, especially some of the
478	! coefficients and indices needed to compute the optical depths
479	! by interpolating data from stored reference atmospheres.
480
481	call setcoef(nlayers, istart, pavel, tavel, tz, tbound, semiss, &
482	coldry, wkl, wbrodl, &
483	laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
484	idrv, dplankbnd_dt, &
485	colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
486	colbrd, fac00, fac01, fac10, fac11, &
487	rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
488	rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
489	rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
490	selffac, selffrac, indself, forfac, forfrac, indfor, &
491	minorfrac, scaleminor, scaleminorn2, indminor)
492
493	! Calculate the gaseous optical depths and Planck fractions for
494	! each longwave spectral band.
495
496	call taumol(nlayers, pavel, wx, coldry, &
497	laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
498	colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
499	colbrd, fac00, fac01, fac10, fac11, &
500	rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
501	rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
502	rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
503	selffac, selffrac, indself, forfac, forfrac, indfor, &
504	minorfrac, scaleminor, scaleminorn2, indminor, &
505	fracs, taug)
506
507	! Combine gaseous and aerosol optical depths, if aerosol active
508	if (iaer .eq. 0) then
509	do k = 1, nlayers
510	do ig = 1, ngptlw
511	taut(k,ig) = taug(k,ig)
512	enddo
513	enddo
514	elseif (iaer .eq. 10) then
515	do k = 1, nlayers
516	do ig = 1, ngptlw
517	taut(k,ig) = taug(k,ig) + taua(k,ngb(ig))
518	enddo
519	enddo
520	endif
521
522	! Call the radiative transfer routine.
523	! Either routine can be called to do clear sky calculation. If clouds
524	! are present, then select routine based on cloud overlap assumption
525	! to be used. Clear sky calculation is done simultaneously.
526
527	if (icld .eq. 1) then
528	call rtrn(nlayers, istart, iend, iout, pz, semiss, ncbands, &
529	cldfrac, taucloud, planklay, planklev, plankbnd, &
530	pwvcm, fracs, taut, &
531	totuflux, totdflux, fnet, htr, &
532	totuclfl, totdclfl, fnetc, htrc, &
533	idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
534	else
535	call rtrnmr(nlayers, istart, iend, iout, pz, semiss, ncbands, &
536	cldfrac, taucloud, planklay, planklev, plankbnd, &
537	pwvcm, fracs, taut, &
538	totuflux, totdflux, fnet, htr, &
539	totuclfl, totdclfl, fnetc, htrc, &
540	idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
541	endif
542
543	! Transfer up and down fluxes and heating rate to output arrays.
544	! Vertical indexing goes from bottom to top; reverse here for GCM if necessary.
545
546	do k = 0, nlayers
547	uflx(iplon,k+1) = totuflux(k)
548	dflx(iplon,k+1) = totdflux(k)
549	uflxc(iplon,k+1) = totuclfl(k)
550	dflxc(iplon,k+1) = totdclfl(k)
551	enddo
552	do k = 0, nlayers-1
553	hr(iplon,k+1) = htr(k)
554	hrc(iplon,k+1) = htrc(k)
555	enddo
556
557	! If idrv=1 option is active, transfer upward flux derivatives to output arrays.
558
559	if (idrv .eq. 1) then
560	do k = 0, nlayers
561	duflx_dt(iplon,k+1) = dtotuflux_dt(k)
562	duflxc_dt(iplon,k+1) = dtotuclfl_dt(k)
563	enddo
564	endif
565
566	! End longitude/column loop
567	enddo
568
569	end subroutine rrtmg_lw
570
571	!***************************************************************************
572	subroutine inatm (iplon, nlay, icld, iaer, &
573	play, plev, tlay, tlev, tsfc, h2ovmr, &
574	o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr, cfc11vmr, cfc12vmr, &
575	cfc22vmr, ccl4vmr, emis, inflglw, iceflglw, liqflglw, &
576	cldfr, taucld, cicewp, cliqwp, reice, reliq, tauaer, &
577	nlayers, pavel, pz, tavel, tz, tbound, semiss, coldry, &
578	wkl, wbrodl, wx, pwvcm, inflag, iceflag, liqflag, &
579	cldfrac, tauc, ciwp, clwp, rei, rel, taua)
580	!***************************************************************************
581	!
582	! Input atmospheric profile from GCM, and prepare it for use in RRTMG_LW.
583	! Set other RRTMG_LW input parameters.
584	!
585	!***************************************************************************
586
587	! --------- Modules ----------
588
589	use parrrtm, only : nbndlw, ngptlw, nmol, maxxsec, mxmol
590	use rrlw_con, only: fluxfac, heatfac, oneminus, pi, grav, avogad
591	use rrlw_wvn, only: ng, nspa, nspb, wavenum1, wavenum2, delwave, ixindx
592
593	! ------- Declarations -------
594
595	! ----- Input -----
596	! Note: All volume mixing ratios are in dimensionless units of mole fraction obtained
597	! by scaling mass mixing ratio (g/g) with the appropriate molecular weights (g/mol)
598	integer(kind=im), intent(in) :: iplon ! column loop index
599	integer(kind=im), intent(in) :: nlay ! Number of model layers
600	integer(kind=im), intent(in) :: icld ! clear/cloud and cloud overlap flag
601	integer(kind=im), intent(in) :: iaer ! aerosol option flag
602
603	real(kind=rb), intent(in) :: play(:,:) ! Layer pressures (hPa, mb)
604	! Dimensions: (ncol,nlay)
605	real(kind=rb), intent(in) :: plev(:,:) ! Interface pressures (hPa, mb)
606	! Dimensions: (ncol,nlay+1)
607	real(kind=rb), intent(in) :: tlay(:,:) ! Layer temperatures (K)
608	! Dimensions: (ncol,nlay)
609	real(kind=rb), intent(in) :: tlev(:,:) ! Interface temperatures (K)
610	! Dimensions: (ncol,nlay+1)
611	real(kind=rb), intent(in) :: tsfc(:) ! Surface temperature (K)
612	! Dimensions: (ncol)
613	real(kind=rb), intent(in) :: h2ovmr(:,:) ! H2O volume mixing ratio
614	! Dimensions: (ncol,nlay)
615	real(kind=rb), intent(in) :: o3vmr(:,:) ! O3 volume mixing ratio
616	! Dimensions: (ncol,nlay)
617	real(kind=rb), intent(in) :: co2vmr(:,:) ! CO2 volume mixing ratio
618	! Dimensions: (ncol,nlay)
619	real(kind=rb), intent(in) :: ch4vmr(:,:) ! Methane volume mixing ratio
620	! Dimensions: (ncol,nlay)
621	real(kind=rb), intent(in) :: n2ovmr(:,:) ! Nitrous oxide volume mixing ratio
622	! Dimensions: (ncol,nlay)
623	real(kind=rb), intent(in) :: o2vmr(:,:) ! Oxygen volume mixing ratio
624	! Dimensions: (ncol,nlay)
625	real(kind=rb), intent(in) :: cfc11vmr(:,:) ! CFC11 volume mixing ratio
626	! Dimensions: (ncol,nlay)
627	real(kind=rb), intent(in) :: cfc12vmr(:,:) ! CFC12 volume mixing ratio
628	! Dimensions: (ncol,nlay)
629	real(kind=rb), intent(in) :: cfc22vmr(:,:) ! CFC22 volume mixing ratio
630	! Dimensions: (ncol,nlay)
631	real(kind=rb), intent(in) :: ccl4vmr(:,:) ! CCL4 volume mixing ratio
632	! Dimensions: (ncol,nlay)
633	real(kind=rb), intent(in) :: emis(:,:) ! Surface emissivity
634	! Dimensions: (ncol,nbndlw)
635
636	integer(kind=im), intent(in) :: inflglw ! Flag for cloud optical properties
637	integer(kind=im), intent(in) :: iceflglw ! Flag for ice particle specification
638	integer(kind=im), intent(in) :: liqflglw ! Flag for liquid droplet specification
639
640	real(kind=rb), intent(in) :: cldfr(:,:) ! Cloud fraction
641	! Dimensions: (ncol,nlay)
642	real(kind=rb), intent(in) :: cicewp(:,:) ! Cloud ice water path (g/m2)
643	! Dimensions: (ncol,nlay)
644	real(kind=rb), intent(in) :: cliqwp(:,:) ! Cloud liquid water path (g/m2)
645	! Dimensions: (ncol,nlay)
646	real(kind=rb), intent(in) :: reice(:,:) ! Cloud ice effective size (microns)
647	! Dimensions: (ncol,nlay)
648	real(kind=rb), intent(in) :: reliq(:,:) ! Cloud water drop effective radius (microns)
649	! Dimensions: (ncol,nlay)
650	real(kind=rb), intent(in) :: taucld(:,:,:) ! In-cloud optical depth
651	! Dimensions: (nbndlw,ncol,nlay)
652	real(kind=rb), intent(in) :: tauaer(:,:,:) ! Aerosol optical depth
653	! Dimensions: (ncol,nlay,nbndlw)
654
655	! ----- Output -----
656	! Atmosphere
657	integer(kind=im), intent(out) :: nlayers ! number of layers
658
659	real(kind=rb), intent(out) :: pavel(:) ! layer pressures (mb)
660	! Dimensions: (nlay)
661	real(kind=rb), intent(out) :: tavel(:) ! layer temperatures (K)
662	! Dimensions: (nlay)
663	real(kind=rb), intent(out) :: pz(0:) ! level (interface) pressures (hPa, mb)
664	! Dimensions: (0:nlay)
665	real(kind=rb), intent(out) :: tz(0:) ! level (interface) temperatures (K)
666	! Dimensions: (0:nlay)
667	real(kind=rb), intent(out) :: tbound ! surface temperature (K)
668	real(kind=rb), intent(out) :: coldry(:) ! dry air column density (mol/cm2)
669	! Dimensions: (nlay)
670	real(kind=rb), intent(out) :: wbrodl(:) ! broadening gas column density (mol/cm2)
671	! Dimensions: (nlay)
672	real(kind=rb), intent(out) :: wkl(:,:) ! molecular amounts (mol/cm-2)
673	! Dimensions: (mxmol,nlay)
674	real(kind=rb), intent(out) :: wx(:,:) ! cross-section amounts (mol/cm-2)
675	! Dimensions: (maxxsec,nlay)
676	real(kind=rb), intent(out) :: pwvcm ! precipitable water vapor (cm)
677	real(kind=rb), intent(out) :: semiss(:) ! lw surface emissivity
678	! Dimensions: (nbndlw)
679
680	! Atmosphere/clouds - cldprop
681	integer(kind=im), intent(out) :: inflag ! flag for cloud property method
682	integer(kind=im), intent(out) :: iceflag ! flag for ice cloud properties
683	integer(kind=im), intent(out) :: liqflag ! flag for liquid cloud properties
684
685	real(kind=rb), intent(out) :: cldfrac(:) ! layer cloud fraction
686	! Dimensions: (nlay)
687	real(kind=rb), intent(out) :: ciwp(:) ! cloud ice water path
688	! Dimensions: (nlay)
689	real(kind=rb), intent(out) :: clwp(:) ! cloud liquid water path
690	! Dimensions: (nlay)
691	real(kind=rb), intent(out) :: rel(:) ! cloud liquid particle effective radius (microns)
692	! Dimensions: (nlay)
693	real(kind=rb), intent(out) :: rei(:) ! cloud ice particle effective size (microns)
694	! Dimensions: (nlay)
695	real(kind=rb), intent(out) :: tauc(:,:) ! in-cloud optical depth
696	! Dimensions: (nbndlw,nlay)
697	real(kind=rb), intent(out) :: taua(:,:) ! aerosol optical depth
698	! Dimensions: (nlay,nbndlw)
699
700
701	! ----- Local -----
702	real(kind=rb), parameter :: amd = 28.9660_rb ! Effective molecular weight of dry air (g/mol)
703	real(kind=rb), parameter :: amw = 18.0160_rb ! Molecular weight of water vapor (g/mol)
704	! real(kind=rb), parameter :: amc = 44.0098_rb ! Molecular weight of carbon dioxide (g/mol)
705	! real(kind=rb), parameter :: amo = 47.9998_rb ! Molecular weight of ozone (g/mol)
706	! real(kind=rb), parameter :: amo2 = 31.9999_rb ! Molecular weight of oxygen (g/mol)
707	! real(kind=rb), parameter :: amch4 = 16.0430_rb ! Molecular weight of methane (g/mol)
708	! real(kind=rb), parameter :: amn2o = 44.0128_rb ! Molecular weight of nitrous oxide (g/mol)
709	! real(kind=rb), parameter :: amc11 = 137.3684_rb ! Molecular weight of CFC11 (g/mol) - CCL3F
710	! real(kind=rb), parameter :: amc12 = 120.9138_rb ! Molecular weight of CFC12 (g/mol) - CCL2F2
711	! real(kind=rb), parameter :: amc22 = 86.4688_rb ! Molecular weight of CFC22 (g/mol) - CHCLF2
712	! real(kind=rb), parameter :: amcl4 = 153.823_rb ! Molecular weight of CCL4 (g/mol) - CCL4
713
714	! Set molecular weight ratios (for converting mmr to vmr)
715	! e.g. h2ovmr = h2ommr * amdw)
716	real(kind=rb), parameter :: amdw = 1.607793_rb ! Molecular weight of dry air / water vapor
717	real(kind=rb), parameter :: amdc = 0.658114_rb ! Molecular weight of dry air / carbon dioxide
718	real(kind=rb), parameter :: amdo = 0.603428_rb ! Molecular weight of dry air / ozone
719	real(kind=rb), parameter :: amdm = 1.805423_rb ! Molecular weight of dry air / methane
720	real(kind=rb), parameter :: amdn = 0.658090_rb ! Molecular weight of dry air / nitrous oxide
721	real(kind=rb), parameter :: amdo2 = 0.905140_rb ! Molecular weight of dry air / oxygen
722	real(kind=rb), parameter :: amdc1 = 0.210852_rb ! Molecular weight of dry air / CFC11
723	real(kind=rb), parameter :: amdc2 = 0.239546_rb ! Molecular weight of dry air / CFC12
724
725	integer(kind=im) :: isp, l, ix, n, imol, ib ! Loop indices
726	real(kind=rb) :: amm, amttl, wvttl, wvsh, summol
727
728	! Add one to nlayers here to include extra model layer at top of atmosphere
729	nlayers = nlay
730
731	! Initialize all molecular amounts and cloud properties to zero here, then pass input amounts
732	! into RRTM arrays below.
733
734	wkl(:,:) = 0.0_rb
735	wx(:,:) = 0.0_rb
736	cldfrac(:) = 0.0_rb
737	tauc(:,:) = 0.0_rb
738	ciwp(:) = 0.0_rb
739	clwp(:) = 0.0_rb
740	rei(:) = 0.0_rb
741	rel(:) = 0.0_rb
742	taua(:,:) = 0.0_rb
743	amttl = 0.0_rb
744	wvttl = 0.0_rb
745
746	! Set surface temperature.
747	tbound = tsfc(iplon)
748
749	! Install input GCM arrays into RRTMG_LW arrays for pressure, temperature,
750	! and molecular amounts.
751	! Pressures are input in mb, or are converted to mb here.
752	! Molecular amounts are input in volume mixing ratio, or are converted from
753	! mass mixing ratio (or specific humidity for h2o) to volume mixing ratio
754	! here. These are then converted to molecular amount (molec/cm2) below.
755	! The dry air column COLDRY (in molec/cm2) is calculated from the level
756	! pressures, pz (in mb), based on the hydrostatic equation and includes a
757	! correction to account for h2o in the layer. The molecular weight of moist
758	! air (amm) is calculated for each layer.
759	! Note: In RRTMG, layer indexing goes from bottom to top, and coding below
760	! assumes GCM input fields are also bottom to top. Input layer indexing
761	! from GCM fields should be reversed here if necessary.
762
763	pz(0) = plev(iplon,1)
764	tz(0) = tlev(iplon,1)
765	do l = 1, nlayers
766	pavel(l) = play(iplon,l)
767	tavel(l) = tlay(iplon,l)
768	pz(l) = plev(iplon,l+1)
769	tz(l) = tlev(iplon,l+1)
770	! For h2o input in vmr:
771	wkl(1,l) = h2ovmr(iplon,l)
772	! For h2o input in mmr:
773	! wkl(1,l) = h2o(iplon,l)*amdw
774	! For h2o input in specific humidity;
775	! wkl(1,l) = (h2o(iplon,l)/(1._rb - h2o(iplon,l)))*amdw
776	wkl(2,l) = co2vmr(iplon,l)
777	wkl(3,l) = o3vmr(iplon,l)
778	wkl(4,l) = n2ovmr(iplon,l)
779	wkl(6,l) = ch4vmr(iplon,l)
780	wkl(7,l) = o2vmr(iplon,l)
781	amm = (1._rb - wkl(1,l)) * amd + wkl(1,l) * amw
782	coldry(l) = (pz(l-1)-pz(l)) * 1.e3_rb * avogad / &
783	(1.e2_rb * grav * amm * (1._rb + wkl(1,l)))
784	enddo
785
786	! Set cross section molecule amounts from input; convert to vmr if necessary
787	do l=1, nlayers
788	wx(1,l) = ccl4vmr(iplon,l)
789	wx(2,l) = cfc11vmr(iplon,l)
790	wx(3,l) = cfc12vmr(iplon,l)
791	wx(4,l) = cfc22vmr(iplon,l)
792	enddo
793
794	! The following section can be used to set values for an additional layer (from
795	! the GCM top level to 1.e-4 mb) for improved calculation of TOA fluxes.
796	! Temperature and molecular amounts in the extra model layer are set to
797	! their values in the top GCM model layer, though these can be modified
798	! here if necessary.
799	! If this feature is utilized, increase nlayers by one above, limit the two
800	! loops above to (nlayers-1), and set the top most (extra) layer values here.
801
802	! pavel(nlayers) = 0.5_rb * pz(nlayers-1)
803	! tavel(nlayers) = tavel(nlayers-1)
804	! pz(nlayers) = 1.e-4_rb
805	! tz(nlayers-1) = 0.5_rb * (tavel(nlayers)+tavel(nlayers-1))
806	! tz(nlayers) = tz(nlayers-1)
807	! wkl(1,nlayers) = wkl(1,nlayers-1)
808	! wkl(2,nlayers) = wkl(2,nlayers-1)
809	! wkl(3,nlayers) = wkl(3,nlayers-1)
810	! wkl(4,nlayers) = wkl(4,nlayers-1)
811	! wkl(6,nlayers) = wkl(6,nlayers-1)
812	! wkl(7,nlayers) = wkl(7,nlayers-1)
813	! amm = (1._rb - wkl(1,nlayers-1)) * amd + wkl(1,nlayers-1) * amw
814	! coldry(nlayers) = (pz(nlayers-1)) * 1.e3_rb * avogad / &
815	! (1.e2_rb * grav * amm * (1._rb + wkl(1,nlayers-1)))
816	! wx(1,nlayers) = wx(1,nlayers-1)
817	! wx(2,nlayers) = wx(2,nlayers-1)
818	! wx(3,nlayers) = wx(3,nlayers-1)
819	! wx(4,nlayers) = wx(4,nlayers-1)
820
821	! At this point all moleculular amounts in wkl and wx are in volume mixing ratio;
822	! convert to molec/cm2 based on coldry for use in rrtm. also, compute precipitable
823	! water vapor for diffusivity angle adjustments in rtrn and rtrnmr.
824
825	do l = 1, nlayers
826	summol = 0.0_rb
827	do imol = 2, nmol
828	summol = summol + wkl(imol,l)
829	enddo
830	wbrodl(l) = coldry(l) * (1._rb - summol)
831	do imol = 1, nmol
832	wkl(imol,l) = coldry(l) * wkl(imol,l)
833	enddo
834	amttl = amttl + coldry(l)+wkl(1,l)
835	wvttl = wvttl + wkl(1,l)
836	do ix = 1,maxxsec
837	if (ixindx(ix) .ne. 0) then
838	wx(ixindx(ix),l) = coldry(l) * wx(ix,l) * 1.e-20_rb
839	endif
840	enddo
841	enddo
842
843	wvsh = (amw * wvttl) / (amd * amttl)
844	pwvcm = wvsh * (1.e3_rb * pz(0)) / (1.e2_rb * grav)
845
846	! Set spectral surface emissivity for each longwave band.
847
848	do n=1,nbndlw
849	semiss(n) = emis(iplon,n)
850	! semiss(n) = 1.0_rb
851	enddo
852
853	! Transfer aerosol optical properties to RRTM variable;
854	! modify to reverse layer indexing here if necessary.
855
856	if (iaer .ge. 1) then
857	do l = 1, nlayers
858	do ib = 1, nbndlw
859	taua(l,ib) = tauaer(iplon,l,ib)
860	enddo
861	enddo
862	endif
863
864	! Transfer cloud fraction and cloud optical properties to RRTM variables,
865	! modify to reverse layer indexing here if necessary.
866
867	if (icld .ge. 1) then
868	inflag = inflglw
869	iceflag = iceflglw
870	liqflag = liqflglw
871
872	! Move incoming GCM cloud arrays to RRTMG cloud arrays.
873	! For GCM input, incoming reice is defined based on selected ice parameterization (inflglw)
874
875	do l = 1, nlayers
876	cldfrac(l) = cldfr(iplon,l)
877	ciwp(l) = cicewp(iplon,l)
878	clwp(l) = cliqwp(iplon,l)
879	rei(l) = reice(iplon,l)
880	rel(l) = reliq(iplon,l)
881	do n=1,nbndlw
882	tauc(n,l) = taucld(n,iplon,l)
883	! ssac(n,l) = ssacld(n,iplon,l)
884	! asmc(n,l) = asmcld(n,iplon,l)
885	enddo
886	enddo
887
888	! If an extra layer is being used in RRTMG, set all cloud properties to zero in the extra layer.
889
890	! cldfrac(nlayers) = 0.0_rb
891	! tauc(:nbndlw,nlayers) = 0.0_rb
892	! ciwp(nlayers) = 0.0_rb
893	! clwp(nlayers) = 0.0_rb
894	! rei(nlayers) = 0.0_rb
895	! rel(nlayers) = 0.0_rb
896	! taua(nlayers,:) = 0.0_rb
897
898	endif
899
900	end subroutine inatm
901
902	end module rrtmg_lw_rad
903

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

Context Navigation

source: palm/trunk/LIB/rrtmg/rrtmg_lw_rad.nomcica.f90 @ 4318

Download in other formats: