Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model.f90 @ 1698

Last change on this file since 1698 was 1692, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 68.2 KB

Line
1	!> @file radiation_model.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2015 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: radiation_model.f90 1692 2015-10-26 16:29:17Z raasch $
26	!
27	! 1691 2015-10-26 16:17:44Z maronga
28	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
29	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
30	! Added output of radiative heating rates.
31	!
32	! 1682 2015-10-07 23:56:08Z knoop
33	! Code annotations made doxygen readable
34	!
35	! 1606 2015-06-29 10:43:37Z maronga
36	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
37	! Note, however, that RRTMG cannot be used without netCDF.
38	!
39	! 1590 2015-05-08 13:56:27Z maronga
40	! Bugfix: definition of character strings requires same length for all elements
41	!
42	! 1587 2015-05-04 14:19:01Z maronga
43	! Added albedo class for snow
44	!
45	! 1585 2015-04-30 07:05:52Z maronga
46	! Added support for RRTMG
47	!
48	! 1571 2015-03-12 16:12:49Z maronga
49	! Added missing KIND attribute. Removed upper-case variable names
50	!
51	! 1551 2015-03-03 14:18:16Z maronga
52	! Added support for data output. Various variables have been renamed. Added
53	! interface for different radiation schemes (currently: clear-sky, constant, and
54	! RRTM (not yet implemented).
55	!
56	! 1496 2014-12-02 17:25:50Z maronga
57	! Initial revision
58	!
59	!
60	! Description:
61	! ------------
62	!> Radiation models and interfaces
63	!> @todo move variable definitions used in init_radiation only to the subroutine
64	!> as they are no longer required after initialization.
65	!> @todo Output of full column vertical profiles used in RRTMG
66	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
67	!> @todo Adapt for use with topography
68	!>
69	!> @note Many variables have a leading dummy dimension (0:0) in order to
70	!> match the assume-size shape expected by the RRTMG model.
71	!------------------------------------------------------------------------------!
72	MODULE radiation_model_mod
73
74
75	USE arrays_3d, &
76	ONLY: dzw, hyp, pt, q, ql, zw
77
78	USE cloud_parameters, &
79	ONLY: cp, l_d_cp, nc_const, rho_l, sigma_gc
80
81	USE constants, &
82	ONLY: pi
83
84	USE control_parameters, &
85	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
86	large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface, &
87	surface_pressure, time_since_reference_point
88
89	USE indices, &
90	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
91
92	USE kinds
93
94	#if defined ( __netcdf )
95	USE netcdf
96	#endif
97
98	USE netcdf_control, &
99	ONLY: dots_label, dots_num, dots_unit
100
101	#if defined ( __rrtmg )
102	USE parrrsw, &
103	ONLY: naerec, nbndsw
104
105	USE parrrtm, &
106	ONLY: nbndlw
107
108	USE rrtmg_lw_init, &
109	ONLY: rrtmg_lw_ini
110
111	USE rrtmg_sw_init, &
112	ONLY: rrtmg_sw_ini
113
114	USE rrtmg_lw_rad, &
115	ONLY: rrtmg_lw
116
117	USE rrtmg_sw_rad, &
118	ONLY: rrtmg_sw
119	#endif
120
121
122
123	IMPLICIT NONE
124
125	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
126
127	!
128	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
129	CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/ &
130	'user defined ', & ! 0
131	'ocean ', & ! 1
132	'mixed farming, tall grassland ', & ! 2
133	'tall/medium grassland ', & ! 3
134	'evergreen shrubland ', & ! 4
135	'short grassland/meadow/shrubland ', & ! 5
136	'evergreen needleleaf forest ', & ! 6
137	'mixed deciduous evergreen forest ', & ! 7
138	'deciduous forest ', & ! 8
139	'tropical evergreen broadleaved forest', & ! 9
140	'medium/tall grassland/woodland ', & ! 10
141	'desert, sandy ', & ! 11
142	'desert, rocky ', & ! 12
143	'tundra ', & ! 13
144	'land ice ', & ! 14
145	'sea ice ', & ! 15
146	'snow ' & ! 16
147	/)
148
149	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
150	day, & !< current day of the year
151	day_init = 172, & !< day of the year at model start (21/06)
152	dots_rad = 0 !< starting index for timeseries output
153
154
155
156
157
158
159	LOGICAL :: constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
160	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
161	lw_radiation = .TRUE., & !< flag parameter indicing whether longwave radiation shall be calculated
162	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
163	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
164	sw_radiation = .TRUE. !< flag parameter indicing whether shortwave radiation shall be calculated
165
166
167	REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp, & !< inverse of seconds per hour (1/3600)
168	d_hours_day = 0.0416666666667_wp, & !< inverse of hours per day (1/24)
169	sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
170	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
171
172	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
173	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
174	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
175	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
176	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
177	decl_1, & !< declination coef. 1
178	decl_2, & !< declination coef. 2
179	decl_3, & !< declination coef. 3
180	dt_radiation = 0.0_wp, & !< radiation model timestep
181	emissivity = 0.98_wp, & !< NAMELIST surface emissivity
182	lambda = 0.0_wp, & !< longitude in degrees
183	lon = 0.0_wp, & !< longitude in radians
184	lat = 0.0_wp, & !< latitude in radians
185	net_radiation = 0.0_wp, & !< net radiation at surface
186	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
187	sky_trans, & !< sky transmissivity
188	time_radiation = 0.0_wp, & !< time since last call of radiation code
189	time_utc, & !< current time in UTC
190	time_utc_init = 43200.0_wp !< UTC time at model start (noon)
191
192	REAL(wp), DIMENSION(0:0) :: zenith !< solar zenith angle
193
194	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
195	alpha, & !< surface broadband albedo (used for clear-sky scheme)
196	rad_net, & !< net radiation at the surface
197	rad_net_av !< average of rad_net
198
199	!
200	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
201	!-- (shortwave, longwave, broadband): sw, lw, bb,
202	REAL(wp), DIMENSION(0:2,1:16), PARAMETER :: albedo_pars = RESHAPE( (/&
203	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
204	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
205	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
206	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
207	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
208	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
209	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
210	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
211	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
212	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
213	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
214	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
215	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
216	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
217	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
218	0.95_wp, 0.70_wp, 0.82_wp & ! 16
219	/), (/ 3, 16 /) )
220
221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
222	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
223	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
224	rad_lw_hr, & !< longwave radiation heating rate (K/s)
225	rad_lw_hr_av, & !< average of rad_sw_hr
226	rad_lw_in, & !< incoming longwave radiation (W/m2)
227	rad_lw_in_av, & !< average of rad_lw_in
228	rad_lw_out, & !< outgoing longwave radiation (W/m2)
229	rad_lw_out_av, & !< average of rad_lw_out
230	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
231	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
232	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
233	rad_sw_hr_av, & !< average of rad_sw_hr
234	rad_sw_in, & !< incoming shortwave radiation (W/m2)
235	rad_sw_in_av, & !< average of rad_sw_in
236	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
237	rad_sw_out_av !< average of rad_sw_out
238
239
240	!
241	!-- Variables and parameters used in RRTMG only
242	#if defined ( __rrtmg )
243	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
244
245
246	!
247	!-- Flag parameters for RRTMGS (should not be changed)
248	INTEGER(iwp), PARAMETER :: rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
249	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
250	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
251	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
252	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
253	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
254
255	!
256	!-- The following variables should be only changed with care, as this will
257	!-- require further setting of some variables, which is currently not
258	!-- implemented (aerosols, ice phase).
259	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
260	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
261	rrtm_iaer = 0, & !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
262	rrtm_idrv = 1 !< longwave upward flux calculation option (0,1)
263
264	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
265
266	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
267
268	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
269	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
270	rrtm_tsfc, & !< dummy array for storing surface temperature
271	t_snd !< actual temperature from sounding data (hPa)
272
273	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
274	aldir, & !< longwave direct albedo solar angle of 60Â°
275	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
276	asdir, & !< shortwave direct albedo solar angle of 60Â°
277	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
278	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
279	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
280	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
281	rrtm_ch4vmr, & !< CH4 volume mixing ratio
282	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
283	rrtm_cldfr, & !< cloud fraction (0,1)
284	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
285	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
286	rrtm_emis, & !< surface emissivity (0-1)
287	rrtm_h2ovmr, & !< H2O volume mixing ratio
288	rrtm_n2ovmr, & !< N2O volume mixing ratio
289	rrtm_o2vmr, & !< O2 volume mixing ratio
290	rrtm_o3vmr, & !< O3 volume mixing ratio
291	rrtm_play, & !< pressure layers (hPa, zu-grid)
292	rrtm_plev, & !< pressure layers (hPa, zw-grid)
293	rrtm_reice, & !< cloud ice effective radius (microns)
294	rrtm_reliq, & !< cloud water drop effective radius (microns)
295	rrtm_tlay, & !< actual temperature (K, zu-grid)
296	rrtm_tlev, & !< actual temperature (K, zw-grid)
297	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
298	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
299	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
300	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
301	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
302	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
303	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
304	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
305	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
306	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
307	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
308	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
309	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
310	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
311
312	!
313	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
314	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
315	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
316	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
317	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
318	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
319	rrtm_aldir, & !< surface albedo for longwave direct radiation
320	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
321	rrtm_asdir, & !< surface albedo for shortwave direct radiation
322	rrtm_lw_tauaer, & !< lw aerosol optical depth
323	rrtm_lw_taucld, & !< lw in-cloud optical depth
324	rrtm_sw_taucld, & !< sw in-cloud optical depth
325	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
326	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
327	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
328	rrtm_sw_tauaer, & !< sw aerosol optical depth
329	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
330	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
331	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
332
333	#endif
334
335	INTERFACE init_radiation
336	MODULE PROCEDURE init_radiation
337	END INTERFACE init_radiation
338
339	INTERFACE radiation_clearsky
340	MODULE PROCEDURE radiation_clearsky
341	END INTERFACE radiation_clearsky
342
343	INTERFACE radiation_rrtmg
344	MODULE PROCEDURE radiation_rrtmg
345	END INTERFACE radiation_rrtmg
346
347	INTERFACE radiation_tendency
348	MODULE PROCEDURE radiation_tendency
349	MODULE PROCEDURE radiation_tendency_ij
350	END INTERFACE radiation_tendency
351
352	SAVE
353
354	PRIVATE
355
356	PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
357	albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad, &
358	dt_radiation, emissivity, force_radiation_call, init_radiation, &
359	lambda, lw_radiation, net_radiation, rad_net, rad_net_av, radiation,&
360	radiation_clearsky, radiation_rrtmg, radiation_scheme, &
361	radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out, &
362	rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, &
363	rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av, &
364	rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, &
365	skip_time_do_radiation, sw_radiation, time_radiation, time_utc_init
366
367
368	#if defined ( __rrtmg )
369	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv, &
370	rrtm_lwuflx_dt
371	#endif
372
373	CONTAINS
374
375	!------------------------------------------------------------------------------!
376	! Description:
377	! ------------
378	!> Initialization of the radiation model
379	!------------------------------------------------------------------------------!
380	SUBROUTINE init_radiation
381
382	IMPLICIT NONE
383
384	!
385	!-- Allocate array for storing the surface net radiation
386	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
387	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
388	rad_net = 0.0_wp
389	ENDIF
390
391	!
392	!-- Fix net radiation in case of radiation_scheme = 'constant'
393	IF ( radiation_scheme == 'constant' ) THEN
394	rad_net = net_radiation
395	radiation = .FALSE.
396	!
397	!-- Calculate orbital constants
398	ELSE
399	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
400	decl_2 = 2.0_wp * pi / 365.0_wp
401	decl_3 = decl_2 * 81.0_wp
402	lat = phi * pi / 180.0_wp
403	lon = lambda * pi / 180.0_wp
404	ENDIF
405
406
407	IF ( radiation_scheme == 'clear-sky' ) THEN
408
409	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
410
411	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
412	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
413	ENDIF
414	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
415	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
416	ENDIF
417
418	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
419	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
420	ENDIF
421	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
422	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
423	ENDIF
424
425	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
426	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
427	ENDIF
428	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
429	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
430	ENDIF
431
432	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
433	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
434	ENDIF
435	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
436	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
437	ENDIF
438
439	rad_sw_in = 0.0_wp
440	rad_sw_out = 0.0_wp
441	rad_lw_in = 0.0_wp
442	rad_lw_out = 0.0_wp
443
444	!
445	!-- Overwrite albedo if manually set in parameter file
446	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
447	albedo = albedo_pars(2,albedo_type)
448	ENDIF
449
450	alpha = albedo
451
452	!
453	!-- Initialization actions for RRTMG
454	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
455	#if defined ( __rrtmg )
456	!
457	!-- Allocate albedos
458	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
459	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
460	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
461	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
462	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
463	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
464	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
465	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
466
467	IF ( albedo_type /= 0 ) THEN
468	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
469	albedo_lw_dif = albedo_pars(0,albedo_type)
470	albedo_lw_dir = albedo_lw_dif
471	ENDIF
472	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
473	albedo_sw_dif = albedo_pars(1,albedo_type)
474	albedo_sw_dir = albedo_sw_dif
475	ENDIF
476	ENDIF
477
478	aldif(:,:) = albedo_lw_dif
479	aldir(:,:) = albedo_lw_dir
480	asdif(:,:) = albedo_sw_dif
481	asdir(:,:) = albedo_sw_dir
482	!
483	!-- Calculate initial values of current (cosine of) the zenith angle and
484	!-- whether the sun is up
485	CALL calc_zenith
486	!
487	!-- Calculate initial surface albedo
488	IF ( .NOT. constant_albedo ) THEN
489	CALL calc_albedo
490	ELSE
491	rrtm_aldif(0,:,:) = aldif(:,:)
492	rrtm_aldir(0,:,:) = aldir(:,:)
493	rrtm_asdif(0,:,:) = asdif(:,:)
494	rrtm_asdir(0,:,:) = asdir(:,:)
495	ENDIF
496
497	!
498	!-- Allocate surface emissivity
499	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
500	rrtm_emis = emissivity
501
502	!
503	!-- Allocate 3d arrays of radiative fluxes and heating rates
504	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
505	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
506	rad_sw_in = 0.0_wp
507	ENDIF
508
509	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
510	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
511	ENDIF
512
513	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
514	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
515	rad_sw_out = 0.0_wp
516	ENDIF
517
518	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
519	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
520	ENDIF
521
522	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
523	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
524	rad_sw_hr = 0.0_wp
525	ENDIF
526
527	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
528	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
529	rad_sw_hr_av = 0.0_wp
530	ENDIF
531
532	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
533	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
534	rad_sw_cs_hr = 0.0_wp
535	ENDIF
536
537	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
538	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
539	rad_sw_cs_hr_av = 0.0_wp
540	ENDIF
541
542	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
543	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
544	rad_lw_in = 0.0_wp
545	ENDIF
546
547	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
548	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
549	ENDIF
550
551	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
552	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
553	rad_lw_out = 0.0_wp
554	ENDIF
555
556	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
557	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
558	ENDIF
559
560	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
561	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
562	rad_lw_hr = 0.0_wp
563	ENDIF
564
565	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
566	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
567	rad_lw_hr_av = 0.0_wp
568	ENDIF
569
570	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
571	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
572	rad_lw_cs_hr = 0.0_wp
573	ENDIF
574
575	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
576	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
577	rad_lw_cs_hr_av = 0.0_wp
578	ENDIF
579
580	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
581	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
582	rad_sw_cs_in = 0.0_wp
583	rad_sw_cs_out = 0.0_wp
584
585	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
586	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
587	rad_lw_cs_in = 0.0_wp
588	rad_lw_cs_out = 0.0_wp
589
590	!
591	!-- Allocate dummy array for storing surface temperature
592	ALLOCATE ( rrtm_tsfc(1) )
593
594	!
595	!-- Initialize RRTMG
596	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
597	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
598
599	!
600	!-- Set input files for RRTMG
601	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
602	IF ( .NOT. snd_exists ) THEN
603	rrtm_input_file = "rrtmg_lw.nc"
604	ENDIF
605
606	!
607	!-- Read vertical layers for RRTMG from sounding data
608	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
609	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
610	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
611	CALL read_sounding_data
612
613	!
614	!-- Read trace gas profiles from file. This routine provides
615	!-- the rrtm_ arrays (1:nzt_rad+1)
616	CALL read_trace_gas_data
617	#endif
618	ENDIF
619
620	!
621	!-- Perform user actions if required
622	CALL user_init_radiation
623
624
625	!
626	!-- Add timeseries for radiation model
627	dots_rad = dots_num + 1
628	dots_num = dots_num + 5
629
630	dots_label(dots_rad+1) = "rad_net"
631	dots_label(dots_rad+2) = "rad_lw_in"
632	dots_label(dots_rad+3) = "rad_lw_out"
633	dots_label(dots_rad+4) = "rad_sw_in"
634	dots_label(dots_rad+5) = "rad_sw_out"
635	dots_unit(dots_rad:dots_rad+4) = "W/m2"
636
637	!
638	!-- Output of albedos is only required for RRTMG
639	IF ( radiation_scheme == 'rrtmg' ) THEN
640	dots_num = dots_num + 4
641	dots_label(dots_rad+5) = "rrtm_aldif"
642	dots_label(dots_rad+6) = "rrtm_aldir"
643	dots_label(dots_rad+7) = "rrtm_asdif"
644	dots_label(dots_rad+8) = "rrtm_asdir"
645	dots_unit(dots_num+5:dots_num+8) = ""
646
647	ENDIF
648
649	!
650	!-- Calculate radiative fluxes at model start
651	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
652	IF ( radiation_scheme == 'clear-sky' ) THEN
653	CALL radiation_clearsky
654	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
655	CALL radiation_rrtmg
656	ENDIF
657	ENDIF
658
659	RETURN
660
661	END SUBROUTINE init_radiation
662
663
664	!------------------------------------------------------------------------------!
665	! Description:
666	! ------------
667	!> A simple clear sky radiation model
668	!------------------------------------------------------------------------------!
669	SUBROUTINE radiation_clearsky
670
671	USE indices, &
672	ONLY: nbgp
673
674	IMPLICIT NONE
675
676	INTEGER(iwp) :: i, j, k !< loop indices
677	REAL(wp) :: exn, & !< Exner functions at surface
678	exn_1, & !< Exner functions at first grid level
679	pt_1 !< potential temperature at first grid level
680
681	!
682	!-- Calculate current zenith angle
683	CALL calc_zenith
684
685	!
686	!-- Calculate sky transmissivity
687	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
688
689	!
690	!-- Calculate value of the Exner function
691	exn = (surface_pressure / 1000.0_wp )**0.286_wp
692	!
693	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
694	!-- point
695	DO i = nxl, nxr
696	DO j = nys, nyn
697	k = nzb_s_inner(j,i)
698
699	exn_1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
700
701	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
702	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
703	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
704
705	IF ( cloud_physics ) THEN
706	pt_1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)
707	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt_1 * exn_1)**4
708	ELSE
709	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn_1)**4
710	ENDIF
711
712	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
713	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
714
715	ENDDO
716	ENDDO
717
718	CALL exchange_horiz_2d( rad_lw_in, nbgp )
719	CALL exchange_horiz_2d( rad_lw_out, nbgp )
720	CALL exchange_horiz_2d( rad_sw_in, nbgp )
721	CALL exchange_horiz_2d( rad_sw_out, nbgp )
722	CALL exchange_horiz_2d( rad_net, nbgp )
723
724	RETURN
725
726	END SUBROUTINE radiation_clearsky
727
728
729	!------------------------------------------------------------------------------!
730	! Description:
731	! ------------
732	!> Implementation of the RRTMG radiation_scheme
733	!------------------------------------------------------------------------------!
734	SUBROUTINE radiation_rrtmg
735
736	USE indices, &
737	ONLY: nbgp
738
739	USE particle_attributes, &
740	ONLY: grid_particles, number_of_particles, particles, &
741	particle_advection_start, prt_count
742
743	IMPLICIT NONE
744
745	#if defined ( __rrtmg )
746
747	INTEGER(iwp) :: i, j, k, n !< loop indices
748
749	REAL(wp) :: s_r2, & !< weighted sum over all droplets with r^2
750	s_r3 !< weighted sum over all droplets with r^3
751
752	!
753	!-- Calculate current (cosine of) zenith angle and whether the sun is up
754	CALL calc_zenith
755	!
756	!-- Calculate surface albedo
757	IF ( .NOT. constant_albedo ) THEN
758	CALL calc_albedo
759	ENDIF
760
761	!
762	!-- Prepare input data for RRTMG
763
764	!
765	!-- In case of large scale forcing with surface data, calculate new pressure
766	!-- profile. nzt_rad might be modified by these calls and all required arrays
767	!-- will then be re-allocated
768	IF ( large_scale_forcing .AND. lsf_surf ) THEN
769	CALL read_sounding_data
770	CALL read_trace_gas_data
771	ENDIF
772	!
773	!-- Loop over all grid points
774	DO i = nxl, nxr
775	DO j = nys, nyn
776
777	!
778	!-- Prepare profiles of temperature and H2O volume mixing ratio
779	rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure &
780	/ 1000.0_wp )**0.286_wp
781
782	DO k = nzb+1, nzt+1
783	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
784	)*0.286_wp + l_d_cp ql(k,j,i)
785	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
786
787	ENDDO
788
789	!
790	!-- Avoid temperature/humidity jumps at the top of the LES domain by
791	!-- linear interpolation from nzt+2 to nzt+7
792	DO k = nzt+2, nzt+7
793	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
794	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
795	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
796	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
797
798	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
799	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
800	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
801	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
802
803	ENDDO
804
805	!-- Linear interpolate to zw grid
806	DO k = nzb+2, nzt+8
807	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
808	rrtm_tlay(0,k-1)) &
809	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
810	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
811	ENDDO
812
813
814	!
815	!-- Calculate liquid water path and cloud fraction for each column.
816	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
817	rrtm_cldfr = 0.0_wp
818	rrtm_reliq = 0.0_wp
819	rrtm_cliqwp = 0.0_wp
820	rrtm_icld = 0
821
822	DO k = nzb+1, nzt+1
823	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
824	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
825	* 100.0_wp / g
826
827	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
828	rrtm_cldfr(0,k) = 1.0_wp
829	IF ( rrtm_icld == 0 ) rrtm_icld = 1
830
831	!
832	!-- Calculate cloud droplet effective radius
833	IF ( cloud_physics ) THEN
834	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
835	* rho_surface &
836	/ ( 4.0_wp * pi * nc_const * rho_l ) &
837	)**0.33333333333333_wp &
838	* EXP( LOG( sigma_gc )**2 )
839
840	ELSEIF ( cloud_droplets ) THEN
841	number_of_particles = prt_count(k,j,i)
842
843	IF (number_of_particles <= 0) CYCLE
844	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
845	s_r2 = 0.0_wp
846	s_r3 = 0.0_wp
847
848	DO n = 1, number_of_particles
849	IF ( particles(n)%particle_mask ) THEN
850	s_r2 = s_r2 + particles(n)%radius*2 &
851	particles(n)%weight_factor
852	s_r3 = s_r3 + particles(n)%radius*3 &
853	particles(n)%weight_factor
854	ENDIF
855	ENDDO
856
857	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
858
859	ENDIF
860
861	!
862	!-- Limit effective radius
863	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
864	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
865	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
866	ENDIF
867	ENDIF
868	ENDDO
869
870	!
871	!-- Set surface temperature
872	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
873
874	IF ( lw_radiation ) THEN
875	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
876	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
877	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
878	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
879	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
880	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
881	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
882	rrtm_reliq , rrtm_lw_tauaer, &
883	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
884	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
885	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
886
887	!
888	!-- Save fluxes
889	DO k = nzb, nzt+1
890	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
891	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
892	ENDDO
893
894	!
895	!-- Save heating rates (convert from K/d to K/h)
896	DO k = nzb+1, nzt+1
897	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k) * d_hours_day
898	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day
899	ENDDO
900
901	ENDIF
902
903	IF ( sw_radiation .AND. sun_up ) THEN
904	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
905	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
906	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
907	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
908	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
909	0.0_wp , day , solar_constant, rrtm_inflgsw,&
910	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
911	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
912	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
913	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
914	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
915	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
916
917	!
918	!-- Save fluxes
919	DO k = nzb, nzt+1
920	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
921	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
922	ENDDO
923
924	!
925	!-- Save heating rates (convert from K/d to K/s)
926	DO k = nzb+1, nzt+1
927	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
928	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
929	ENDDO
930
931	ENDIF
932
933	!
934	!-- Calculate surface net radiation
935	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
936	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
937
938	ENDDO
939	ENDDO
940
941	CALL exchange_horiz( rad_lw_in, nbgp )
942	CALL exchange_horiz( rad_lw_out, nbgp )
943	CALL exchange_horiz( rad_lw_hr, nbgp )
944	CALL exchange_horiz( rad_lw_cs_hr, nbgp )
945
946	CALL exchange_horiz( rad_sw_in, nbgp )
947	CALL exchange_horiz( rad_sw_out, nbgp )
948	CALL exchange_horiz( rad_sw_hr, nbgp )
949	CALL exchange_horiz( rad_sw_cs_hr, nbgp )
950
951	CALL exchange_horiz_2d( rad_net, nbgp )
952	#endif
953
954	END SUBROUTINE radiation_rrtmg
955
956
957	!------------------------------------------------------------------------------!
958	! Description:
959	! ------------
960	!> Calculate the cosine of the zenith angle (variable is called zenith)
961	!------------------------------------------------------------------------------!
962	SUBROUTINE calc_zenith
963
964	IMPLICIT NONE
965
966	REAL(wp) :: declination, & !< solar declination angle
967	hour_angle !< solar hour angle
968	!
969	!-- Calculate current day and time based on the initial values and simulation
970	!-- time
971	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
972	/ 86400.0_wp ), KIND=iwp)
973	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
974
975
976	!
977	!-- Calculate solar declination and hour angle
978	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
979	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
980
981	!
982	!-- Calculate zenith angle
983	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
984	* COS(hour_angle)
985	zenith(0) = MAX(0.0_wp,zenith(0))
986
987	!
988	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
989	IF ( zenith(0) > 0.0_wp ) THEN
990	sun_up = .TRUE.
991	ELSE
992	sun_up = .FALSE.
993	END IF
994
995	END SUBROUTINE calc_zenith
996
997	#if defined ( __rrtmg ) && defined ( __netcdf )
998	!------------------------------------------------------------------------------!
999	! Description:
1000	! ------------
1001	!> Calculates surface albedo components based on Briegleb (1992) and
1002	!> Briegleb et al. (1986)
1003	!------------------------------------------------------------------------------!
1004	SUBROUTINE calc_albedo
1005
1006	IMPLICIT NONE
1007
1008	IF ( sun_up ) THEN
1009	!
1010	!-- Ocean
1011	IF ( albedo_type == 1 ) THEN
1012	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
1013	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
1014	* ( zenith(0) - 0.5_wp ) &
1015	* ( zenith(0) - 1.0_wp )
1016	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1017	!
1018	!-- Snow
1019	ELSEIF ( albedo_type == 16 ) THEN
1020	IF ( zenith(0) < 0.5_wp ) THEN
1021	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
1022	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1023	* zenith(0))) - 1.0_wp
1024	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
1025	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1026	* zenith(0))) - 1.0_wp
1027
1028	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1029	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1030	ELSE
1031	rrtm_aldir(0,:,:) = aldif
1032	rrtm_asdir(0,:,:) = asdif
1033	ENDIF
1034	!
1035	!-- Sea ice
1036	ELSEIF ( albedo_type == 15 ) THEN
1037	rrtm_aldir(0,:,:) = aldif
1038	rrtm_asdir(0,:,:) = asdif
1039	!
1040	!-- Land surfaces
1041	ELSE
1042	SELECT CASE ( albedo_type )
1043
1044	!
1045	!-- Surface types with strong zenith dependence
1046	CASE ( 1, 2, 3, 4, 11, 12, 13 )
1047	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
1048	(1.0_wp + 0.8_wp * zenith(0))
1049	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
1050	(1.0_wp + 0.8_wp * zenith(0))
1051	!
1052	!-- Surface types with weak zenith dependence
1053	CASE ( 5, 6, 7, 8, 9, 10, 14 )
1054	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
1055	(1.0_wp + 0.2_wp * zenith(0))
1056	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
1057	(1.0_wp + 0.2_wp * zenith(0))
1058
1059	CASE DEFAULT
1060
1061	END SELECT
1062	ENDIF
1063	!
1064	!-- Diffusive albedo is taken from Table 2
1065	rrtm_aldif(0,:,:) = aldif
1066	rrtm_asdif(0,:,:) = asdif
1067
1068	ELSE
1069
1070	rrtm_aldir(0,:,:) = 0.0_wp
1071	rrtm_asdir(0,:,:) = 0.0_wp
1072	rrtm_aldif(0,:,:) = 0.0_wp
1073	rrtm_asdif(0,:,:) = 0.0_wp
1074	ENDIF
1075	END SUBROUTINE calc_albedo
1076
1077	!------------------------------------------------------------------------------!
1078	! Description:
1079	! ------------
1080	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1081	!------------------------------------------------------------------------------!
1082	SUBROUTINE read_sounding_data
1083
1084	USE netcdf_control
1085
1086	IMPLICIT NONE
1087
1088	INTEGER(iwp) :: id, & !< NetCDF id of input file
1089	id_dim_zrad, & !< pressure level id in the NetCDF file
1090	id_var, & !< NetCDF variable id
1091	k, & !< loop index
1092	nz_snd, & !< number of vertical levels in the sounding data
1093	nz_snd_start, & !< start vertical index for sounding data to be used
1094	nz_snd_end !< end vertical index for souding data to be used
1095
1096	REAL(wp) :: t_surface !< actual surface temperature
1097
1098	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1099	t_snd_tmp !< temporary temperature profile (sounding)
1100
1101	!
1102	!-- In case of updates, deallocate arrays first (sufficient to check one
1103	!-- array as the others are automatically allocated). This is required
1104	!-- because nzt_rad might change during the update
1105	IF ( ALLOCATED ( hyp_snd ) ) THEN
1106	DEALLOCATE( hyp_snd )
1107	DEALLOCATE( t_snd )
1108	DEALLOCATE( q_snd )
1109	DEALLOCATE ( rrtm_play )
1110	DEALLOCATE ( rrtm_plev )
1111	DEALLOCATE ( rrtm_tlay )
1112	DEALLOCATE ( rrtm_tlev )
1113
1114	DEALLOCATE ( rrtm_h2ovmr )
1115	DEALLOCATE ( rrtm_cicewp )
1116	DEALLOCATE ( rrtm_cldfr )
1117	DEALLOCATE ( rrtm_cliqwp )
1118	DEALLOCATE ( rrtm_reice )
1119	DEALLOCATE ( rrtm_reliq )
1120	DEALLOCATE ( rrtm_lw_taucld )
1121	DEALLOCATE ( rrtm_lw_tauaer )
1122
1123	DEALLOCATE ( rrtm_lwdflx )
1124	DEALLOCATE ( rrtm_lwdflxc )
1125	DEALLOCATE ( rrtm_lwuflx )
1126	DEALLOCATE ( rrtm_lwuflxc )
1127	DEALLOCATE ( rrtm_lwuflx_dt )
1128	DEALLOCATE ( rrtm_lwuflxc_dt )
1129	DEALLOCATE ( rrtm_lwhr )
1130	DEALLOCATE ( rrtm_lwhrc )
1131
1132	DEALLOCATE ( rrtm_sw_taucld )
1133	DEALLOCATE ( rrtm_sw_ssacld )
1134	DEALLOCATE ( rrtm_sw_asmcld )
1135	DEALLOCATE ( rrtm_sw_fsfcld )
1136	DEALLOCATE ( rrtm_sw_tauaer )
1137	DEALLOCATE ( rrtm_sw_ssaaer )
1138	DEALLOCATE ( rrtm_sw_asmaer )
1139	DEALLOCATE ( rrtm_sw_ecaer )
1140
1141	DEALLOCATE ( rrtm_swdflx )
1142	DEALLOCATE ( rrtm_swdflxc )
1143	DEALLOCATE ( rrtm_swuflx )
1144	DEALLOCATE ( rrtm_swuflxc )
1145	DEALLOCATE ( rrtm_swhr )
1146	DEALLOCATE ( rrtm_swhrc )
1147
1148	ENDIF
1149
1150	!
1151	!-- Open file for reading
1152	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1153	CALL handle_netcdf_error( 'netcdf', 549 )
1154
1155	!
1156	!-- Inquire dimension of z axis and save in nz_snd
1157	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1158	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1159	CALL handle_netcdf_error( 'netcdf', 551 )
1160
1161	!
1162	! !-- Allocate temporary array for storing pressure data
1163	ALLOCATE( hyp_snd_tmp(nzb+1:nz_snd) )
1164	hyp_snd_tmp = 0.0_wp
1165
1166
1167	!-- Read pressure from file
1168	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1169	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
1170	count = (/nz_snd/) )
1171	CALL handle_netcdf_error( 'netcdf', 552 )
1172
1173	!
1174	!-- Allocate temporary array for storing temperature data
1175	ALLOCATE( t_snd_tmp(nzb+1:nz_snd) )
1176	t_snd_tmp = 0.0_wp
1177
1178	!
1179	!-- Read temperature from file
1180	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1181	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
1182	count = (/nz_snd/) )
1183	CALL handle_netcdf_error( 'netcdf', 553 )
1184
1185	!
1186	!-- Calculate start of sounding data
1187	nz_snd_start = nz_snd + 1
1188	nz_snd_end = nz_snd_end
1189
1190	!
1191	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
1192	!-- in Pa, hyp_snd in hPa).
1193	DO k = 1, nz_snd
1194	IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
1195	nz_snd_start = k
1196	EXIT
1197	END IF
1198	END DO
1199
1200	IF ( nz_snd_start <= nz_snd ) THEN
1201	nz_snd_end = nz_snd - 1
1202	END IF
1203
1204
1205	!
1206	!-- Calculate of total grid points for RRTMG calculations
1207	nzt_rad = nzt + nz_snd_end - nz_snd_start + 2
1208
1209	!
1210	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
1211	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1212	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1213	ALLOCATE( t_snd(nzb+1:nzt_rad) )
1214	ALLOCATE( q_snd(nzb+1:nzt_rad) )
1215	hyp_snd = 0.0_wp
1216	t_snd = 0.0_wp
1217	q_snd = 0.0_wp
1218
1219	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start:nz_snd_end)
1220	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start:nz_snd_end)
1221
1222	DEALLOCATE ( hyp_snd_tmp )
1223	DEALLOCATE ( t_snd_tmp )
1224
1225	nc_stat = NF90_CLOSE( id )
1226
1227	!
1228	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
1229	!-- top of the LES domain. This routine does not consider horizontal or
1230	!-- vertical variability of pressure and temperature
1231	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
1232	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
1233
1234	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1235	DO k = nzb+1, nzt+1
1236	rrtm_play(0,k) = hyp(k) * 0.01_wp
1237	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1238	t_surface )**(1.0_wp/0.286_wp)
1239	ENDDO
1240
1241	DO k = nzt+2, nzt_rad
1242	rrtm_play(0,k) = hyp_snd(k)
1243	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1244	ENDDO
1245	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
1246	1.5 * hyp_snd(nzt_rad) &
1247	- 0.5 * hyp_snd(nzt_rad-1) )
1248	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
1249	0.25_wp * rrtm_plev(0,nzt_rad+1) )
1250
1251	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1252
1253	!
1254	!-- Calculate temperature/humidity levels at top of the LES domain.
1255	!-- Currently, the temperature is taken from sounding data (might lead to a
1256	!-- temperature jump at interface. To do: Humidity is currently not
1257	!-- calculated above the LES domain.
1258	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
1259	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
1260	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1261
1262	DO k = nzt+8, nzt_rad
1263	rrtm_tlay(0,k) = t_snd(k)
1264	rrtm_h2ovmr(0,k) = q_snd(k)
1265	ENDDO
1266	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
1267	- rrtm_tlay(0,nzt_rad-1)
1268	DO k = nzt+9, nzt_rad+1
1269	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
1270	- rrtm_tlay(0,k-1)) &
1271	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
1272	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1273	ENDDO
1274	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1275
1276	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
1277	- rrtm_tlev(0,nzt_rad)
1278	!
1279	!-- Allocate remaining RRTMG arrays
1280	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1281	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1282	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1283	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1284	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1285	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1286	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1287	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1288	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1289	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1290	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1291	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1292	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1293	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1294	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
1295
1296	!
1297	!-- The ice phase is currently not considered in PALM
1298	rrtm_cicewp = 0.0_wp
1299	rrtm_reice = 0.0_wp
1300
1301	!
1302	!-- Set other parameters (move to NAMELIST parameters in the future)
1303	rrtm_lw_tauaer = 0.0_wp
1304	rrtm_lw_taucld = 0.0_wp
1305	rrtm_sw_taucld = 0.0_wp
1306	rrtm_sw_ssacld = 0.0_wp
1307	rrtm_sw_asmcld = 0.0_wp
1308	rrtm_sw_fsfcld = 0.0_wp
1309	rrtm_sw_tauaer = 0.0_wp
1310	rrtm_sw_ssaaer = 0.0_wp
1311	rrtm_sw_asmaer = 0.0_wp
1312	rrtm_sw_ecaer = 0.0_wp
1313
1314
1315	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
1316	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
1317	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
1318	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1319	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1320	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1321
1322	rrtm_swdflx = 0.0_wp
1323	rrtm_swuflx = 0.0_wp
1324	rrtm_swhr = 0.0_wp
1325	rrtm_swuflxc = 0.0_wp
1326	rrtm_swdflxc = 0.0_wp
1327	rrtm_swhrc = 0.0_wp
1328
1329	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
1330	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
1331	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
1332	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1333	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1334	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1335
1336	rrtm_lwdflx = 0.0_wp
1337	rrtm_lwuflx = 0.0_wp
1338	rrtm_lwhr = 0.0_wp
1339	rrtm_lwuflxc = 0.0_wp
1340	rrtm_lwdflxc = 0.0_wp
1341	rrtm_lwhrc = 0.0_wp
1342
1343	ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
1344	ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
1345
1346	rrtm_lwuflxc_dt = 0.0_wp
1347	rrtm_lwuflxc_dt = 0.0_wp
1348
1349	END SUBROUTINE read_sounding_data
1350
1351
1352	!------------------------------------------------------------------------------!
1353	! Description:
1354	! ------------
1355	!> Read trace gas data from file
1356	!------------------------------------------------------------------------------!
1357	SUBROUTINE read_trace_gas_data
1358
1359	USE netcdf_control
1360	USE rrsw_ncpar
1361
1362	IMPLICIT NONE
1363
1364	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
1365
1366	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: & !< trace gas names
1367	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
1368	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
1369
1370	INTEGER(iwp) :: id, & !< NetCDF id
1371	k, & !< loop index
1372	m, & !< loop index
1373	n, & !< loop index
1374	nabs, & !< number of absorbers
1375	np, & !< number of pressure levels
1376	id_abs, & !< NetCDF id of the respective absorber
1377	id_dim, & !< NetCDF id of asborber's dimension
1378	id_var !< NetCDf id ot the absorber
1379
1380	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
1381
1382
1383	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
1384	rrtm_play_tmp, & !< temporary array for pressure zu-levels
1385	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
1386	trace_path_tmp !< temporary array for storing trace gas path data
1387
1388	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
1389	trace_mls_path, & !< array for storing trace gas path data
1390	trace_mls_tmp !< temporary array for storing trace gas data
1391
1392
1393	!
1394	!-- In case of updates, deallocate arrays first (sufficient to check one
1395	!-- array as the others are automatically allocated)
1396	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
1397	DEALLOCATE ( rrtm_o3vmr )
1398	DEALLOCATE ( rrtm_co2vmr )
1399	DEALLOCATE ( rrtm_ch4vmr )
1400	DEALLOCATE ( rrtm_n2ovmr )
1401	DEALLOCATE ( rrtm_o2vmr )
1402	DEALLOCATE ( rrtm_cfc11vmr )
1403	DEALLOCATE ( rrtm_cfc12vmr )
1404	DEALLOCATE ( rrtm_cfc22vmr )
1405	DEALLOCATE ( rrtm_ccl4vmr )
1406	ENDIF
1407
1408	!
1409	!-- Allocate trace gas profiles
1410	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
1411	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
1412	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
1413	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
1414	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
1415	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
1416	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
1417	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
1418	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
1419
1420	!
1421	!-- Open file for reading
1422	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1423	CALL handle_netcdf_error( 'netcdf', 549 )
1424	!
1425	!-- Inquire dimension ids and dimensions
1426	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
1427	CALL handle_netcdf_error( 'netcdf', 550 )
1428	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
1429	CALL handle_netcdf_error( 'netcdf', 550 )
1430
1431	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
1432	CALL handle_netcdf_error( 'netcdf', 550 )
1433	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
1434	CALL handle_netcdf_error( 'netcdf', 550 )
1435
1436
1437	!
1438	!-- Allocate pressure, and trace gas arrays
1439	ALLOCATE( p_mls(1:np) )
1440	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
1441	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
1442
1443
1444	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1445	CALL handle_netcdf_error( 'netcdf', 550 )
1446	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
1447	CALL handle_netcdf_error( 'netcdf', 550 )
1448
1449	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
1450	CALL handle_netcdf_error( 'netcdf', 550 )
1451	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
1452	CALL handle_netcdf_error( 'netcdf', 550 )
1453
1454
1455	!
1456	!-- Write absorber amounts (mls) to trace_mls
1457	DO n = 1, num_trace_gases
1458	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
1459
1460	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
1461
1462	!
1463	!-- Replace missing values by zero
1464	WHERE ( trace_mls(n,:) > 2.0_wp )
1465	trace_mls(n,:) = 0.0_wp
1466	END WHERE
1467	END DO
1468
1469	DEALLOCATE ( trace_mls_tmp )
1470
1471	nc_stat = NF90_CLOSE( id )
1472	CALL handle_netcdf_error( 'netcdf', 551 )
1473
1474	!
1475	!-- Add extra pressure level for calculations of the trace gas paths
1476	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
1477	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
1478
1479	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
1480	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
1481	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
1482	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
1483	* rrtm_plev(0,nzt_rad+1) )
1484
1485	!
1486	!-- Calculate trace gas path (zero at surface) with interpolation to the
1487	!-- sounding levels
1488	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
1489
1490	trace_mls_path(nzb+1,:) = 0.0_wp
1491
1492	DO k = nzb+2, nzt_rad+2
1493	DO m = 1, num_trace_gases
1494	trace_mls_path(k,m) = trace_mls_path(k-1,m)
1495
1496	!
1497	!-- When the pressure level is higher than the trace gas pressure
1498	!-- level, assume that
1499	IF ( rrtm_plev_tmp(k-1) > p_mls(1) ) THEN
1500
1501	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
1502	* ( rrtm_plev_tmp(k-1) &
1503	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
1504	) / g
1505	ENDIF
1506
1507	!
1508	!-- Integrate for each sounding level from the contributing p_mls
1509	!-- levels
1510	DO n = 2, np
1511	!
1512	!-- Limit p_mls so that it is within the model level
1513	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
1514	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
1515	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
1516	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
1517
1518	IF ( p_mls_l > p_mls_u ) THEN
1519
1520	!
1521	!-- Calculate weights for interpolation
1522	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
1523	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
1524	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
1525
1526	!
1527	!-- Add level to trace gas path
1528	trace_mls_path(k,m) = trace_mls_path(k,m) &
1529	+ ( p_wgt_u * trace_mls(m,n) &
1530	+ p_wgt_l * trace_mls(m,n-1) ) &
1531	* (p_mls_l - p_mls_u) / g
1532	ENDIF
1533	ENDDO
1534
1535	IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN
1536	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
1537	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
1538	- rrtm_plev_tmp(k) &
1539	) / g
1540	ENDIF
1541	ENDDO
1542	ENDDO
1543
1544
1545	!
1546	!-- Prepare trace gas path profiles
1547	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
1548
1549	DO m = 1, num_trace_gases
1550
1551	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
1552	- trace_mls_path(1:nzt_rad+1,m) ) * g &
1553	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
1554	- rrtm_plev_tmp(2:nzt_rad+2) )
1555
1556	!
1557	!-- Save trace gas paths to the respective arrays
1558	SELECT CASE ( TRIM( trace_names(m) ) )
1559
1560	CASE ( 'O3' )
1561
1562	rrtm_o3vmr(0,:) = trace_path_tmp(:)
1563
1564	CASE ( 'CO2' )
1565
1566	rrtm_co2vmr(0,:) = trace_path_tmp(:)
1567
1568	CASE ( 'CH4' )
1569
1570	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
1571
1572	CASE ( 'N2O' )
1573
1574	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
1575
1576	CASE ( 'O2' )
1577
1578	rrtm_o2vmr(0,:) = trace_path_tmp(:)
1579
1580	CASE ( 'CFC11' )
1581
1582	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
1583
1584	CASE ( 'CFC12' )
1585
1586	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
1587
1588	CASE ( 'CFC22' )
1589
1590	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
1591
1592	CASE ( 'CCL4' )
1593
1594	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
1595
1596	CASE DEFAULT
1597
1598	END SELECT
1599
1600	ENDDO
1601
1602	DEALLOCATE ( trace_path_tmp )
1603	DEALLOCATE ( trace_mls_path )
1604	DEALLOCATE ( rrtm_play_tmp )
1605	DEALLOCATE ( rrtm_plev_tmp )
1606	DEALLOCATE ( trace_mls )
1607	DEALLOCATE ( p_mls )
1608
1609	END SUBROUTINE read_trace_gas_data
1610
1611	#endif
1612
1613
1614	!------------------------------------------------------------------------------!
1615	! Description:
1616	! ------------
1617	!> Calculate temperature tendency due to radiative cooling/heating.
1618	!> Cache-optimized version.
1619	!------------------------------------------------------------------------------!
1620	SUBROUTINE radiation_tendency_ij ( i, j, tend )
1621
1622	USE cloud_parameters, &
1623	ONLY: pt_d_t
1624
1625	IMPLICIT NONE
1626
1627	INTEGER(iwp) :: i, j, k !< loop indices
1628
1629	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
1630
1631	#if defined ( __rrtmg )
1632	!
1633	!-- Calculate tendency based on heating rate
1634	DO k = nzb+1, nzt+1
1635	tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) &
1636	* pt_d_t(k) * d_seconds_hour
1637	ENDDO
1638
1639	#endif
1640
1641	END SUBROUTINE radiation_tendency_ij
1642
1643
1644	!------------------------------------------------------------------------------!
1645	! Description:
1646	! ------------
1647	!> Calculate temperature tendency due to radiative cooling/heating.
1648	!> Vector-optimized version
1649	!------------------------------------------------------------------------------!
1650	SUBROUTINE radiation_tendency ( tend )
1651
1652	USE cloud_parameters, &
1653	ONLY: pt_d_t
1654
1655	USE indices, &
1656	ONLY: nxl, nxr, nyn, nys
1657
1658	IMPLICIT NONE
1659
1660	INTEGER(iwp) :: i, j, k !< loop indices
1661
1662	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
1663
1664	#if defined ( __rrtmg )
1665	!
1666	!-- Calculate tendency based on heating rate
1667	DO i = nxl, nxr
1668	DO j = nys, nyn
1669	DO k = nzb+1, nzt+1
1670	tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) &
1671	+ rad_sw_hr(k,j,i) ) * pt_d_t(k) &
1672	* d_seconds_hour
1673	ENDDO
1674	ENDDO
1675	ENDDO
1676	#endif
1677
1678	END SUBROUTINE radiation_tendency
1679
1680	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |