Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 2238

Last change on this file since 2238 was 2233, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 143.6 KB

Line
1	!> @file radiation_model_mod.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
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: radiation_model_mod.f90 2233 2017-05-30 18:08:54Z suehring $
27	!
28	! 2232 2017-05-30 17:47:52Z suehring
29	! Adjustments to new topography concept
30	! Bugfix in read restart
31	!
32	! 2200 2017-04-11 11:37:51Z suehring
33	! Bugfix in call of exchange_horiz_2d and read restart data
34	!
35	! 2163 2017-03-01 13:23:15Z schwenkel
36	! Bugfix in radiation_check_data_output
37	!
38	! 2157 2017-02-22 15:10:35Z suehring
39	! Bugfix in read_restart data
40	!
41	! 2011 2016-09-19 17:29:57Z kanani
42	! Removed CALL of auxiliary SUBROUTINE get_usm_info,
43	! flag urban_surface is now defined in module control_parameters.
44	!
45	! 2007 2016-08-24 15:47:17Z kanani
46	! Added calculation of solar directional vector for new urban surface
47	! model,
48	! accounted for urban_surface model in radiation_check_parameters,
49	! correction of comments for zenith angle.
50	!
51	! 2000 2016-08-20 18:09:15Z knoop
52	! Forced header and separation lines into 80 columns
53	!
54	! 1976 2016-07-27 13:28:04Z maronga
55	! Output of 2D/3D/masked data is now directly done within this module. The
56	! radiation schemes have been simplified for better usability so that
57	! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
58	! the radiation code used.
59	!
60	! 1856 2016-04-13 12:56:17Z maronga
61	! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
62	!
63	! 1853 2016-04-11 09:00:35Z maronga
64	! Added routine for radiation_scheme = constant.
65	!
66	! 1849 2016-04-08 11:33:18Z hoffmann
67	! Adapted for modularization of microphysics
68	!
69	! 1826 2016-04-07 12:01:39Z maronga
70	! Further modularization.
71	!
72	! 1788 2016-03-10 11:01:04Z maronga
73	! Added new albedo class for pavements / roads.
74	!
75	! 1783 2016-03-06 18:36:17Z raasch
76	! palm-netcdf-module removed in order to avoid a circular module dependency,
77	! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
78	! added
79	!
80	! 1757 2016-02-22 15:49:32Z maronga
81	! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
82	! profiles for pressure and temperature above the LES domain.
83	!
84	! 1709 2015-11-04 14:47:01Z maronga
85	! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
86	! corrections
87	!
88	! 1701 2015-11-02 07:43:04Z maronga
89	! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
90	!
91	! 1691 2015-10-26 16:17:44Z maronga
92	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
93	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
94	! Added output of radiative heating rates.
95	!
96	! 1682 2015-10-07 23:56:08Z knoop
97	! Code annotations made doxygen readable
98	!
99	! 1606 2015-06-29 10:43:37Z maronga
100	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
101	! Note, however, that RRTMG cannot be used without netCDF.
102	!
103	! 1590 2015-05-08 13:56:27Z maronga
104	! Bugfix: definition of character strings requires same length for all elements
105	!
106	! 1587 2015-05-04 14:19:01Z maronga
107	! Added albedo class for snow
108	!
109	! 1585 2015-04-30 07:05:52Z maronga
110	! Added support for RRTMG
111	!
112	! 1571 2015-03-12 16:12:49Z maronga
113	! Added missing KIND attribute. Removed upper-case variable names
114	!
115	! 1551 2015-03-03 14:18:16Z maronga
116	! Added support for data output. Various variables have been renamed. Added
117	! interface for different radiation schemes (currently: clear-sky, constant, and
118	! RRTM (not yet implemented).
119	!
120	! 1496 2014-12-02 17:25:50Z maronga
121	! Initial revision
122	!
123	!
124	! Description:
125	! ------------
126	!> Radiation models and interfaces
127	!> @todo move variable definitions used in radiation_init only to the subroutine
128	!> as they are no longer required after initialization.
129	!> @todo Output of full column vertical profiles used in RRTMG
130	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
131	!> @todo Adapt for use with topography
132	!>
133	!> @note Many variables have a leading dummy dimension (0:0) in order to
134	!> match the assume-size shape expected by the RRTMG model.
135	!------------------------------------------------------------------------------!
136	MODULE radiation_model_mod
137
138	USE arrays_3d, &
139	ONLY: dzw, hyp, pt, q, ql, zu, zw
140
141	USE cloud_parameters, &
142	ONLY: cp, l_d_cp, rho_l
143
144	USE constants, &
145	ONLY: pi
146
147	USE control_parameters, &
148	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
149	large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface, &
150	surface_pressure, time_since_reference_point
151
152	USE indices, &
153	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
154	wall_flags_0
155
156	USE kinds
157
158	USE microphysics_mod, &
159	ONLY: nc_const, sigma_gc
160
161	#if defined ( __netcdf )
162	USE NETCDF
163	#endif
164
165	#if defined ( __rrtmg )
166	USE parrrsw, &
167	ONLY: naerec, nbndsw
168
169	USE parrrtm, &
170	ONLY: nbndlw
171
172	USE rrtmg_lw_init, &
173	ONLY: rrtmg_lw_ini
174
175	USE rrtmg_sw_init, &
176	ONLY: rrtmg_sw_ini
177
178	USE rrtmg_lw_rad, &
179	ONLY: rrtmg_lw
180
181	USE rrtmg_sw_rad, &
182	ONLY: rrtmg_sw
183	#endif
184
185
186
187	IMPLICIT NONE
188
189	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
190
191	!
192	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
193	CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/ &
194	'user defined ', & ! 0
195	'ocean ', & ! 1
196	'mixed farming, tall grassland ', & ! 2
197	'tall/medium grassland ', & ! 3
198	'evergreen shrubland ', & ! 4
199	'short grassland/meadow/shrubland ', & ! 5
200	'evergreen needleleaf forest ', & ! 6
201	'mixed deciduous evergreen forest ', & ! 7
202	'deciduous forest ', & ! 8
203	'tropical evergreen broadleaved forest', & ! 9
204	'medium/tall grassland/woodland ', & ! 10
205	'desert, sandy ', & ! 11
206	'desert, rocky ', & ! 12
207	'tundra ', & ! 13
208	'land ice ', & ! 14
209	'sea ice ', & ! 15
210	'snow ', & ! 16
211	'pavement/roads ' & ! 17
212	/)
213
214	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
215	day, & !< current day of the year
216	day_init = 172, & !< day of the year at model start (21/06)
217	dots_rad = 0 !< starting index for timeseries output
218
219	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
220	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
221	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
222	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
223	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
224	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
225	sw_radiation = .TRUE., & !< flag parameter indicing whether shortwave radiation shall be calculated
226	sun_direction = .FALSE. !< flag parameter indicing whether solar direction shall be calculated
227
228
229	REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp, & !< inverse of seconds per hour (1/3600)
230	d_hours_day = 0.0416666666667_wp, & !< inverse of hours per day (1/24)
231	sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
232	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
233
234	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
235	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
236	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
237	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
238	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
239	decl_1, & !< declination coef. 1
240	decl_2, & !< declination coef. 2
241	decl_3, & !< declination coef. 3
242	dt_radiation = 0.0_wp, & !< radiation model timestep
243	emissivity = 0.98_wp, & !< NAMELIST surface emissivity
244	lambda = 0.0_wp, & !< longitude in degrees
245	lon = 0.0_wp, & !< longitude in radians
246	lat = 0.0_wp, & !< latitude in radians
247	net_radiation = 0.0_wp, & !< net radiation at surface
248	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
249	sky_trans, & !< sky transmissivity
250	time_radiation = 0.0_wp, & !< time since last call of radiation code
251	time_utc, & !< current time in UTC
252	time_utc_init = 43200.0_wp !< UTC time at model start (noon)
253
254	REAL(wp), DIMENSION(0:0) :: zenith, & !< cosine of solar zenith angle
255	sun_dir_lat, & !< solar directional vector in latitudes
256	sun_dir_lon !< solar directional vector in longitudes
257
258	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
259	alpha, & !< surface broadband albedo (used for clear-sky scheme)
260	rad_lw_out_change_0, & !< change in LW out due to change in surface temperature
261	rad_net, & !< net radiation at the surface
262	rad_net_av !< average of rad_net
263
264	!
265	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
266	!-- (shortwave, longwave, broadband): sw, lw, bb,
267	REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
268	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
269	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
270	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
271	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
272	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
273	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
274	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
275	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
276	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
277	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
278	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
279	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
280	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
281	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
282	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
283	0.95_wp, 0.70_wp, 0.82_wp, & ! 16
284	0.08_wp, 0.08_wp, 0.08_wp & ! 17
285	/), (/ 3, 17 /) )
286
287	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
288	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
289	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
290	rad_lw_hr, & !< longwave radiation heating rate (K/s)
291	rad_lw_hr_av, & !< average of rad_sw_hr
292	rad_lw_in, & !< incoming longwave radiation (W/m2)
293	rad_lw_in_av, & !< average of rad_lw_in
294	rad_lw_out, & !< outgoing longwave radiation (W/m2)
295	rad_lw_out_av, & !< average of rad_lw_out
296	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
297	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
298	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
299	rad_sw_hr_av, & !< average of rad_sw_hr
300	rad_sw_in, & !< incoming shortwave radiation (W/m2)
301	rad_sw_in_av, & !< average of rad_sw_in
302	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
303	rad_sw_out_av !< average of rad_sw_out
304
305
306	!
307	!-- Variables and parameters used in RRTMG only
308	#if defined ( __rrtmg )
309	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
310
311
312	!
313	!-- Flag parameters for RRTMGS (should not be changed)
314	INTEGER(iwp), PARAMETER :: rrtm_idrv = 1, & !< flag for longwave upward flux calculation option (0,1)
315	rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
316	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
317	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
318	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
319	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
320	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
321
322	!
323	!-- The following variables should be only changed with care, as this will
324	!-- require further setting of some variables, which is currently not
325	!-- implemented (aerosols, ice phase).
326	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
327	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
328	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)
329
330	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
331
332	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
333
334	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
335
336	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
337	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
338	rrtm_tsfc, & !< dummy array for storing surface temperature
339	t_snd !< actual temperature from sounding data (hPa)
340
341	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
342	aldir, & !< longwave direct albedo solar angle of 60Â°
343	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
344	asdir, & !< shortwave direct albedo solar angle of 60Â°
345	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
346	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
347	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
348	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
349	rrtm_ch4vmr, & !< CH4 volume mixing ratio
350	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
351	rrtm_cldfr, & !< cloud fraction (0,1)
352	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
353	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
354	rrtm_emis, & !< surface emissivity (0-1)
355	rrtm_h2ovmr, & !< H2O volume mixing ratio
356	rrtm_n2ovmr, & !< N2O volume mixing ratio
357	rrtm_o2vmr, & !< O2 volume mixing ratio
358	rrtm_o3vmr, & !< O3 volume mixing ratio
359	rrtm_play, & !< pressure layers (hPa, zu-grid)
360	rrtm_plev, & !< pressure layers (hPa, zw-grid)
361	rrtm_reice, & !< cloud ice effective radius (microns)
362	rrtm_reliq, & !< cloud water drop effective radius (microns)
363	rrtm_tlay, & !< actual temperature (K, zu-grid)
364	rrtm_tlev, & !< actual temperature (K, zw-grid)
365	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
366	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
367	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
368	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
369	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
370	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
371	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
372	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
373	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
374	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
375	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
376	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
377	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
378	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
379
380	!
381	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
382	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
383	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
384	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
385	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
386	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
387	rrtm_aldir, & !< surface albedo for longwave direct radiation
388	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
389	rrtm_asdir, & !< surface albedo for shortwave direct radiation
390	rrtm_lw_tauaer, & !< lw aerosol optical depth
391	rrtm_lw_taucld, & !< lw in-cloud optical depth
392	rrtm_sw_taucld, & !< sw in-cloud optical depth
393	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
394	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
395	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
396	rrtm_sw_tauaer, & !< sw aerosol optical depth
397	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
398	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
399	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
400
401	#endif
402
403	INTERFACE radiation_check_data_output
404	MODULE PROCEDURE radiation_check_data_output
405	END INTERFACE radiation_check_data_output
406
407	INTERFACE radiation_check_data_output_pr
408	MODULE PROCEDURE radiation_check_data_output_pr
409	END INTERFACE radiation_check_data_output_pr
410
411	INTERFACE radiation_check_parameters
412	MODULE PROCEDURE radiation_check_parameters
413	END INTERFACE radiation_check_parameters
414
415	INTERFACE radiation_clearsky
416	MODULE PROCEDURE radiation_clearsky
417	END INTERFACE radiation_clearsky
418
419	INTERFACE radiation_constant
420	MODULE PROCEDURE radiation_constant
421	END INTERFACE radiation_constant
422
423	INTERFACE radiation_control
424	MODULE PROCEDURE radiation_control
425	END INTERFACE radiation_control
426
427	INTERFACE radiation_3d_data_averaging
428	MODULE PROCEDURE radiation_3d_data_averaging
429	END INTERFACE radiation_3d_data_averaging
430
431	INTERFACE radiation_data_output_2d
432	MODULE PROCEDURE radiation_data_output_2d
433	END INTERFACE radiation_data_output_2d
434
435	INTERFACE radiation_data_output_3d
436	MODULE PROCEDURE radiation_data_output_3d
437	END INTERFACE radiation_data_output_3d
438
439	INTERFACE radiation_data_output_mask
440	MODULE PROCEDURE radiation_data_output_mask
441	END INTERFACE radiation_data_output_mask
442
443	INTERFACE radiation_define_netcdf_grid
444	MODULE PROCEDURE radiation_define_netcdf_grid
445	END INTERFACE radiation_define_netcdf_grid
446
447	INTERFACE radiation_header
448	MODULE PROCEDURE radiation_header
449	END INTERFACE radiation_header
450
451	INTERFACE radiation_init
452	MODULE PROCEDURE radiation_init
453	END INTERFACE radiation_init
454
455	INTERFACE radiation_parin
456	MODULE PROCEDURE radiation_parin
457	END INTERFACE radiation_parin
458
459	INTERFACE radiation_rrtmg
460	MODULE PROCEDURE radiation_rrtmg
461	END INTERFACE radiation_rrtmg
462
463	INTERFACE radiation_tendency
464	MODULE PROCEDURE radiation_tendency
465	MODULE PROCEDURE radiation_tendency_ij
466	END INTERFACE radiation_tendency
467
468	INTERFACE radiation_read_restart_data
469	MODULE PROCEDURE radiation_read_restart_data
470	END INTERFACE radiation_read_restart_data
471
472	INTERFACE radiation_last_actions
473	MODULE PROCEDURE radiation_last_actions
474	END INTERFACE radiation_last_actions
475
476	SAVE
477
478	PRIVATE
479
480	!
481	!-- Public functions / NEEDS SORTING
482	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
483	radiation_check_parameters, radiation_control, &
484	radiation_header, radiation_init, radiation_parin, &
485	radiation_3d_data_averaging, radiation_tendency, &
486	radiation_data_output_2d, radiation_data_output_3d, &
487	radiation_define_netcdf_grid, radiation_last_actions, &
488	radiation_read_restart_data, radiation_data_output_mask
489
490	!
491	!-- Public variables and constants / NEEDS SORTING
492	PUBLIC dots_rad, dt_radiation, force_radiation_call, &
493	rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in, &
494	rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, &
495	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
496	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
497	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, &
498	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
499	zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, &
500	day_init, time_utc_init
501
502
503	#if defined ( __rrtmg )
504	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
505	#endif
506
507	CONTAINS
508
509
510	!------------------------------------------------------------------------------!
511	! Description:
512	! ------------
513	!> This subroutine controls the calls of the radiation schemes
514	!------------------------------------------------------------------------------!
515	SUBROUTINE radiation_control
516
517
518	IMPLICIT NONE
519
520
521	SELECT CASE ( TRIM( radiation_scheme ) )
522
523	CASE ( 'constant' )
524	CALL radiation_constant
525
526	CASE ( 'clear-sky' )
527	CALL radiation_clearsky
528
529	CASE ( 'rrtmg' )
530	CALL radiation_rrtmg
531
532	CASE DEFAULT
533
534	END SELECT
535
536
537	END SUBROUTINE radiation_control
538
539	!------------------------------------------------------------------------------!
540	! Description:
541	! ------------
542	!> Check data output for radiation model
543	!------------------------------------------------------------------------------!
544	SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
545
546
547	USE control_parameters, &
548	ONLY: data_output, message_string
549
550	IMPLICIT NONE
551
552	CHARACTER (LEN=*) :: unit !<
553	CHARACTER (LEN=*) :: var !<
554
555	INTEGER(iwp) :: i
556	INTEGER(iwp) :: ilen
557	INTEGER(iwp) :: k
558
559	SELECT CASE ( TRIM( var ) )
560
561	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
562	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
563	message_string = '"output of "' // TRIM( var ) // '" requi' // &
564	'res radiation = .TRUE. and ' // &
565	'radiation_scheme = "rrtmg"'
566	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
567	ENDIF
568	unit = 'K/h'
569
570	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
571	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
572	message_string = '"output of "' // TRIM( var ) // '" requi' // &
573	'res radiation = .TRUE. and ' // &
574	'radiation_scheme = "rrtmg"'
575	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
576	ENDIF
577	unit = 'W/m2'
578
579	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
580	'rrtm_asdir*' )
581	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
582	message_string = 'illegal value for data_output: "' // &
583	TRIM( var ) // '" & only 2d-horizontal ' // &
584	'cross sections are allowed for this value'
585	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
586	ENDIF
587	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
588	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
589	TRIM( var ) == 'rrtm_aldir*' .OR. &
590	TRIM( var ) == 'rrtm_asdif*' .OR. &
591	TRIM( var ) == 'rrtm_asdir*' ) &
592	THEN
593	message_string = 'output of "' // TRIM( var ) // '" require'&
594	// 's radiation = .TRUE. and radiation_sch'&
595	// 'eme = "rrtmg"'
596	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
597	ENDIF
598	ENDIF
599
600	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
601	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
602	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
603	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
604	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
605
606	CASE DEFAULT
607	unit = 'illegal'
608
609	END SELECT
610
611
612	END SUBROUTINE radiation_check_data_output
613
614	!------------------------------------------------------------------------------!
615	! Description:
616	! ------------
617	!> Check data output of profiles for radiation model
618	!------------------------------------------------------------------------------!
619	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
620
621	USE arrays_3d, &
622	ONLY: zu
623
624	USE control_parameters, &
625	ONLY: data_output_pr, message_string
626
627	USE indices
628
629	USE profil_parameter
630
631	USE statistics
632
633	IMPLICIT NONE
634
635	CHARACTER (LEN=*) :: unit !<
636	CHARACTER (LEN=*) :: variable !<
637	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
638
639	INTEGER(iwp) :: user_pr_index !<
640	INTEGER(iwp) :: var_count !<
641
642	SELECT CASE ( TRIM( variable ) )
643
644	CASE ( 'rad_net' )
645	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
646	THEN
647	message_string = 'data_output_pr = ' // &
648	TRIM( data_output_pr(var_count) ) // ' is' // &
649	'not available for radiation = .FALSE. or ' //&
650	'radiation_scheme = "constant"'
651	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
652	ELSE
653	dopr_index(var_count) = 101
654	dopr_unit = 'W/m2'
655	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
656	unit = dopr_unit
657	ENDIF
658
659	CASE ( 'rad_lw_in' )
660	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
661	THEN
662	message_string = 'data_output_pr = ' // &
663	TRIM( data_output_pr(var_count) ) // ' is' // &
664	'not available for radiation = .FALSE. or ' //&
665	'radiation_scheme = "constant"'
666	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
667	ELSE
668	dopr_index(var_count) = 102
669	dopr_unit = 'W/m2'
670	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
671	unit = dopr_unit
672	ENDIF
673
674	CASE ( 'rad_lw_out' )
675	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
676	THEN
677	message_string = 'data_output_pr = ' // &
678	TRIM( data_output_pr(var_count) ) // ' is' // &
679	'not available for radiation = .FALSE. or ' //&
680	'radiation_scheme = "constant"'
681	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
682	ELSE
683	dopr_index(var_count) = 103
684	dopr_unit = 'W/m2'
685	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
686	unit = dopr_unit
687	ENDIF
688
689	CASE ( 'rad_sw_in' )
690	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
691	THEN
692	message_string = 'data_output_pr = ' // &
693	TRIM( data_output_pr(var_count) ) // ' is' // &
694	'not available for radiation = .FALSE. or ' //&
695	'radiation_scheme = "constant"'
696	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
697	ELSE
698	dopr_index(var_count) = 104
699	dopr_unit = 'W/m2'
700	hom(:,2,104,:) = SPREAD( zw, 2, statistic_regions+1 )
701	unit = dopr_unit
702	ENDIF
703
704	CASE ( 'rad_sw_out')
705	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
706	THEN
707	message_string = 'data_output_pr = ' // &
708	TRIM( data_output_pr(var_count) ) // ' is' // &
709	'not available for radiation = .FALSE. or ' //&
710	'radiation_scheme = "constant"'
711	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
712	ELSE
713	dopr_index(var_count) = 105
714	dopr_unit = 'W/m2'
715	hom(:,2,105,:) = SPREAD( zw, 2, statistic_regions+1 )
716	unit = dopr_unit
717	ENDIF
718
719	CASE ( 'rad_lw_cs_hr' )
720	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
721	THEN
722	message_string = 'data_output_pr = ' // &
723	TRIM( data_output_pr(var_count) ) // ' is' // &
724	'not available for radiation = .FALSE. or ' //&
725	'radiation_scheme /= "rrtmg"'
726	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
727	ELSE
728	dopr_index(var_count) = 106
729	dopr_unit = 'K/h'
730	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
731	unit = dopr_unit
732	ENDIF
733
734	CASE ( 'rad_lw_hr' )
735	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
736	THEN
737	message_string = 'data_output_pr = ' // &
738	TRIM( data_output_pr(var_count) ) // ' is' // &
739	'not available for radiation = .FALSE. or ' //&
740	'radiation_scheme /= "rrtmg"'
741	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
742	ELSE
743	dopr_index(var_count) = 107
744	dopr_unit = 'K/h'
745	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
746	unit = dopr_unit
747	ENDIF
748
749	CASE ( 'rad_sw_cs_hr' )
750	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
751	THEN
752	message_string = 'data_output_pr = ' // &
753	TRIM( data_output_pr(var_count) ) // ' is' // &
754	'not available for radiation = .FALSE. or ' //&
755	'radiation_scheme /= "rrtmg"'
756	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
757	ELSE
758	dopr_index(var_count) = 108
759	dopr_unit = 'K/h'
760	hom(:,2,108,:) = SPREAD( zu, 2, statistic_regions+1 )
761	unit = dopr_unit
762	ENDIF
763
764	CASE ( 'rad_sw_hr' )
765	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
766	THEN
767	message_string = 'data_output_pr = ' // &
768	TRIM( data_output_pr(var_count) ) // ' is' // &
769	'not available for radiation = .FALSE. or ' //&
770	'radiation_scheme /= "rrtmg"'
771	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
772	ELSE
773	dopr_index(var_count) = 109
774	dopr_unit = 'K/h'
775	hom(:,2,109,:) = SPREAD( zu, 2, statistic_regions+1 )
776	unit = dopr_unit
777	ENDIF
778
779
780	CASE DEFAULT
781	unit = 'illegal'
782
783	END SELECT
784
785
786	END SUBROUTINE radiation_check_data_output_pr
787
788
789	!------------------------------------------------------------------------------!
790	! Description:
791	! ------------
792	!> Check parameters routine for radiation model
793	!------------------------------------------------------------------------------!
794	SUBROUTINE radiation_check_parameters
795
796	USE control_parameters, &
797	ONLY: message_string, topography, urban_surface
798
799
800	IMPLICIT NONE
801
802
803	IF ( radiation_scheme /= 'constant' .AND. &
804	radiation_scheme /= 'clear-sky' .AND. &
805	radiation_scheme /= 'rrtmg' ) THEN
806	message_string = 'unknown radiation_scheme = '// &
807	TRIM( radiation_scheme )
808	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
809	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
810	#if ! defined ( __rrtmg )
811	message_string = 'radiation_scheme = "rrtmg" requires ' // &
812	'compilation of PALM with pre-processor ' // &
813	'directive -D__rrtmg'
814	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
815	#endif
816	#if defined ( __rrtmg ) && ! defined( __netcdf )
817	message_string = 'radiation_scheme = "rrtmg" requires ' // &
818	'the use of NetCDF (preprocessor directive ' // &
819	'-D__netcdf'
820	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
821	#endif
822
823	ENDIF
824
825	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
826	radiation_scheme == 'clear-sky') THEN
827	message_string = 'radiation_scheme = "clear-sky" in combination' // &
828	'with albedo_type = 0 requires setting of albedo'// &
829	' /= 9999999.9'
830	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
831	ENDIF
832
833	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
834	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
835	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
836	) ) THEN
837	message_string = 'radiation_scheme = "rrtmg" in combination' // &
838	'with albedo_type = 0 requires setting of ' // &
839	'albedo_lw_dif /= 9999999.9' // &
840	'albedo_lw_dir /= 9999999.9' // &
841	'albedo_sw_dif /= 9999999.9 and' // &
842	'albedo_sw_dir /= 9999999.9'
843	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
844	ENDIF
845
846	!
847	!-- The following paramter check is temporarily extended by the urban_surface
848	!-- flag, until a better solution comes up to omit this check in case of
849	!-- urban surface model is used.
850	IF ( topography /= 'flat' .AND. .NOT. urban_surface ) THEN
851	message_string = 'radiation scheme cannot be used ' // &
852	'in combination with topography /= "flat"'
853	CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
854	ENDIF
855
856	END SUBROUTINE radiation_check_parameters
857
858
859	!------------------------------------------------------------------------------!
860	! Description:
861	! ------------
862	!> Initialization of the radiation model
863	!------------------------------------------------------------------------------!
864	SUBROUTINE radiation_init
865
866	IMPLICIT NONE
867
868	!
869	!-- Allocate array for storing the surface net radiation
870	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
871	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
872	rad_net = 0.0_wp
873	ENDIF
874
875	!
876	!-- Allocate array for storing the surface net radiation
877	IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) ) THEN
878	ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
879	rad_lw_out_change_0 = 0.0_wp
880	ENDIF
881
882	!
883	!-- Fix net radiation in case of radiation_scheme = 'constant'
884	IF ( radiation_scheme == 'constant' ) THEN
885	rad_net = net_radiation
886	! radiation = .FALSE.
887	!
888	!-- Calculate orbital constants
889	ELSE
890	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
891	decl_2 = 2.0_wp * pi / 365.0_wp
892	decl_3 = decl_2 * 81.0_wp
893	lat = phi * pi / 180.0_wp
894	lon = lambda * pi / 180.0_wp
895	ENDIF
896
897	IF ( radiation_scheme == 'clear-sky' .OR. &
898	radiation_scheme == 'constant') THEN
899
900	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
901
902	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
903	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
904	ENDIF
905	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
906	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
907	ENDIF
908
909	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
910	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
911	ENDIF
912	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
913	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
914	ENDIF
915
916	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
917	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
918	ENDIF
919	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
920	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
921	ENDIF
922
923	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
924	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
925	ENDIF
926	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
927	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
928	ENDIF
929
930	rad_sw_in = 0.0_wp
931	rad_sw_out = 0.0_wp
932	rad_lw_in = 0.0_wp
933	rad_lw_out = 0.0_wp
934
935	!
936	!-- Overwrite albedo if manually set in parameter file
937	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
938	albedo = albedo_pars(2,albedo_type)
939	ENDIF
940
941	alpha = albedo
942
943	!
944	!-- Initialization actions for RRTMG
945	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
946	#if defined ( __rrtmg )
947	!
948	!-- Allocate albedos
949	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
950	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
951	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
952	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
953	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
954	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
955	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
956	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
957
958	IF ( albedo_type /= 0 ) THEN
959	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
960	albedo_lw_dif = albedo_pars(0,albedo_type)
961	albedo_lw_dir = albedo_lw_dif
962	ENDIF
963	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
964	albedo_sw_dif = albedo_pars(1,albedo_type)
965	albedo_sw_dir = albedo_sw_dif
966	ENDIF
967	ENDIF
968
969	aldif(:,:) = albedo_lw_dif
970	aldir(:,:) = albedo_lw_dir
971	asdif(:,:) = albedo_sw_dif
972	asdir(:,:) = albedo_sw_dir
973	!
974	!-- Calculate initial values of current (cosine of) the zenith angle and
975	!-- whether the sun is up
976	CALL calc_zenith
977	!
978	!-- Calculate initial surface albedo
979	IF ( .NOT. constant_albedo ) THEN
980	CALL calc_albedo
981	ELSE
982	rrtm_aldif(0,:,:) = aldif(:,:)
983	rrtm_aldir(0,:,:) = aldir(:,:)
984	rrtm_asdif(0,:,:) = asdif(:,:)
985	rrtm_asdir(0,:,:) = asdir(:,:)
986	ENDIF
987
988	!
989	!-- Allocate surface emissivity
990	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
991	rrtm_emis = emissivity
992
993	!
994	!-- Allocate 3d arrays of radiative fluxes and heating rates
995	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
996	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
997	rad_sw_in = 0.0_wp
998	ENDIF
999
1000	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
1001	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1002	ENDIF
1003
1004	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
1005	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1006	rad_sw_out = 0.0_wp
1007	ENDIF
1008
1009	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
1010	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1011	ENDIF
1012
1013	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
1014	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1015	rad_sw_hr = 0.0_wp
1016	ENDIF
1017
1018	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
1019	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1020	rad_sw_hr_av = 0.0_wp
1021	ENDIF
1022
1023	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
1024	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1025	rad_sw_cs_hr = 0.0_wp
1026	ENDIF
1027
1028	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
1029	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1030	rad_sw_cs_hr_av = 0.0_wp
1031	ENDIF
1032
1033	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
1034	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1035	rad_lw_in = 0.0_wp
1036	ENDIF
1037
1038	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
1039	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1040	ENDIF
1041
1042	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
1043	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1044	rad_lw_out = 0.0_wp
1045	ENDIF
1046
1047	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
1048	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1049	ENDIF
1050
1051	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
1052	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1053	rad_lw_hr = 0.0_wp
1054	ENDIF
1055
1056	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
1057	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1058	rad_lw_hr_av = 0.0_wp
1059	ENDIF
1060
1061	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
1062	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1063	rad_lw_cs_hr = 0.0_wp
1064	ENDIF
1065
1066	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
1067	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1068	rad_lw_cs_hr_av = 0.0_wp
1069	ENDIF
1070
1071	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1072	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1073	rad_sw_cs_in = 0.0_wp
1074	rad_sw_cs_out = 0.0_wp
1075
1076	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1077	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1078	rad_lw_cs_in = 0.0_wp
1079	rad_lw_cs_out = 0.0_wp
1080
1081	!
1082	!-- Allocate dummy array for storing surface temperature
1083	ALLOCATE ( rrtm_tsfc(1) )
1084
1085	!
1086	!-- Initialize RRTMG
1087	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
1088	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
1089
1090	!
1091	!-- Set input files for RRTMG
1092	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
1093	IF ( .NOT. snd_exists ) THEN
1094	rrtm_input_file = "rrtmg_lw.nc"
1095	ENDIF
1096
1097	!
1098	!-- Read vertical layers for RRTMG from sounding data
1099	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
1100	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
1101	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
1102	CALL read_sounding_data
1103
1104	!
1105	!-- Read trace gas profiles from file. This routine provides
1106	!-- the rrtm_ arrays (1:nzt_rad+1)
1107	CALL read_trace_gas_data
1108	#endif
1109	ENDIF
1110
1111	!
1112	!-- Perform user actions if required
1113	CALL user_init_radiation
1114
1115	!
1116	!-- Calculate radiative fluxes at model start
1117	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1118
1119	SELECT CASE ( radiation_scheme )
1120	CASE ( 'rrtmg' )
1121	CALL radiation_rrtmg
1122	CASE ( 'clear-sky' )
1123	CALL radiation_clearsky
1124	CASE ( 'constant' )
1125	CALL radiation_constant
1126	CASE DEFAULT
1127	END SELECT
1128
1129	ENDIF
1130
1131	RETURN
1132
1133	END SUBROUTINE radiation_init
1134
1135
1136	!------------------------------------------------------------------------------!
1137	! Description:
1138	! ------------
1139	!> A simple clear sky radiation model
1140	!------------------------------------------------------------------------------!
1141	SUBROUTINE radiation_clearsky
1142
1143
1144	IMPLICIT NONE
1145
1146	INTEGER(iwp) :: i, j, k !< loop indices
1147	REAL(wp) :: exn, & !< Exner functions at surface
1148	exn1, & !< Exner functions at first grid level
1149	pt1 !< potential temperature at first grid level
1150
1151	!
1152	!-- Calculate current zenith angle
1153	CALL calc_zenith
1154
1155	!
1156	!-- Calculate sky transmissivity
1157	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1158
1159	!
1160	!-- Calculate value of the Exner function
1161	exn = (surface_pressure / 1000.0_wp )**0.286_wp
1162	!
1163	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
1164	!-- point
1165	DO i = nxlg, nxrg
1166	DO j = nysg, nyng
1167	!
1168	!-- Obtain vertical index of topography top
1169	k = MAXLOC( &
1170	MERGE( 1, 0, &
1171	BTEST( wall_flags_0(:,j,i), 12 ) &
1172	), DIM = 1 &
1173	) - 1
1174
1175	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1176
1177	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
1178	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1179	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1180
1181	IF ( cloud_physics ) THEN
1182	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1183	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1184	ELSE
1185	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1186	ENDIF
1187
1188	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
1189	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1190
1191
1192	rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity &
1193	* (pt(k,j,i) * exn) ** 3
1194
1195	ENDDO
1196	ENDDO
1197
1198	END SUBROUTINE radiation_clearsky
1199
1200
1201	!------------------------------------------------------------------------------!
1202	! Description:
1203	! ------------
1204	!> This scheme keeps the prescribed net radiation constant during the run
1205	!------------------------------------------------------------------------------!
1206	SUBROUTINE radiation_constant
1207
1208
1209	IMPLICIT NONE
1210
1211	INTEGER(iwp) :: i, j, k !< loop indices
1212	REAL(wp) :: exn, & !< Exner functions at surface
1213	exn1, & !< Exner functions at first grid level
1214	pt1 !< potential temperature at first grid level
1215
1216	!
1217	!-- Calculate value of the Exner function
1218	exn = (surface_pressure / 1000.0_wp )**0.286_wp
1219	!
1220	!-- Prescribe net radiation and estimate the remaining radiative fluxes
1221	DO i = nxlg, nxrg
1222	DO j = nysg, nyng
1223	!
1224	!-- Obtain vertical index of topography top. So far it is identical to nzb.
1225	k = MAXLOC( &
1226	MERGE( 1, 0, &
1227	BTEST( wall_flags_0(:,j,i), 12 ) &
1228	), DIM = 1 &
1229	) - 1
1230
1231	rad_net(j,i) = net_radiation
1232
1233	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1234
1235	IF ( cloud_physics ) THEN
1236	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1237	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1238	ELSE
1239	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1240	ENDIF
1241
1242	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1243
1244	rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i) &
1245	+ rad_lw_out(0,j,i) ) &
1246	/ ( 1.0_wp - alpha(j,i) )
1247
1248	ENDDO
1249	ENDDO
1250
1251	END SUBROUTINE radiation_constant
1252
1253	!------------------------------------------------------------------------------!
1254	! Description:
1255	! ------------
1256	!> Header output for radiation model
1257	!------------------------------------------------------------------------------!
1258	SUBROUTINE radiation_header ( io )
1259
1260
1261	IMPLICIT NONE
1262
1263	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1264
1265
1266
1267	!
1268	!-- Write radiation model header
1269	WRITE( io, 3 )
1270
1271	IF ( radiation_scheme == "constant" ) THEN
1272	WRITE( io, 4 ) net_radiation
1273	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
1274	WRITE( io, 5 )
1275	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
1276	WRITE( io, 6 )
1277	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
1278	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
1279	ENDIF
1280
1281	IF ( albedo_type == 0 ) THEN
1282	WRITE( io, 7 ) albedo
1283	ELSE
1284	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1285	ENDIF
1286	IF ( constant_albedo ) THEN
1287	WRITE( io, 9 )
1288	ENDIF
1289
1290	IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN
1291	WRITE ( io, 1 ) lambda
1292	WRITE ( io, 2 ) day_init, time_utc_init
1293	ENDIF
1294
1295	WRITE( io, 12 ) dt_radiation
1296
1297
1298	1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr')
1299	2 FORMAT (' Day of the year at model start : day_init = ',I3 &
1300	/' UTC time at model start : time_utc_init = ',F7.1' s')
1301	3 FORMAT (//' Radiation model information:'/ &
1302	' ----------------------------'/)
1303	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
1304	// 'W/m**2')
1305	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', &
1306	' default)')
1307	6 FORMAT (' --> RRTMG scheme is used')
1308	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
1309	8 FORMAT (/' Albedo is set for land surface type: ', A)
1310	9 FORMAT (/' --> Albedo is fixed during the run')
1311	10 FORMAT (/' --> Longwave radiation is disabled')
1312	11 FORMAT (/' --> Shortwave radiation is disabled.')
1313	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
1314
1315
1316	END SUBROUTINE radiation_header
1317
1318
1319	!------------------------------------------------------------------------------!
1320	! Description:
1321	! ------------
1322	!> Parin for &radiation_par for radiation model
1323	!------------------------------------------------------------------------------!
1324	SUBROUTINE radiation_parin
1325
1326
1327	IMPLICIT NONE
1328
1329	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
1330
1331	NAMELIST /radiation_par/ albedo, albedo_type, albedo_lw_dir, &
1332	albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1333	constant_albedo, day_init, dt_radiation, &
1334	lambda, lw_radiation, net_radiation, &
1335	radiation_scheme, skip_time_do_radiation, &
1336	sw_radiation, time_utc_init, &
1337	unscheduled_radiation_calls
1338
1339	line = ' '
1340
1341	!
1342	!-- Try to find radiation model package
1343	REWIND ( 11 )
1344	line = ' '
1345	DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1346	READ ( 11, '(A)', END=10 ) line
1347	ENDDO
1348	BACKSPACE ( 11 )
1349
1350	!
1351	!-- Read user-defined namelist
1352	READ ( 11, radiation_par )
1353
1354	!
1355	!-- Set flag that indicates that the radiation model is switched on
1356	radiation = .TRUE.
1357
1358	10 CONTINUE
1359
1360
1361	END SUBROUTINE radiation_parin
1362
1363
1364	!------------------------------------------------------------------------------!
1365	! Description:
1366	! ------------
1367	!> Implementation of the RRTMG radiation_scheme
1368	!------------------------------------------------------------------------------!
1369	SUBROUTINE radiation_rrtmg
1370
1371	USE indices, &
1372	ONLY: nbgp
1373
1374	USE particle_attributes, &
1375	ONLY: grid_particles, number_of_particles, particles, &
1376	particle_advection_start, prt_count
1377
1378	IMPLICIT NONE
1379
1380	#if defined ( __rrtmg )
1381
1382	INTEGER(iwp) :: i, j, k, n !< loop indices
1383
1384	REAL(wp) :: s_r2, & !< weighted sum over all droplets with r^2
1385	s_r3 !< weighted sum over all droplets with r^3
1386
1387	!
1388	!-- Calculate current (cosine of) zenith angle and whether the sun is up
1389	CALL calc_zenith
1390	!
1391	!-- Calculate surface albedo
1392	IF ( .NOT. constant_albedo ) THEN
1393	CALL calc_albedo
1394	ENDIF
1395
1396	!
1397	!-- Prepare input data for RRTMG
1398
1399	!
1400	!-- In case of large scale forcing with surface data, calculate new pressure
1401	!-- profile. nzt_rad might be modified by these calls and all required arrays
1402	!-- will then be re-allocated
1403	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1404	CALL read_sounding_data
1405	CALL read_trace_gas_data
1406	ENDIF
1407	!
1408	!-- Loop over all grid points
1409	DO i = nxl, nxr
1410	DO j = nys, nyn
1411
1412	!
1413	!-- Prepare profiles of temperature and H2O volume mixing ratio
1414	rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure &
1415	/ 1000.0_wp )**0.286_wp
1416
1417	DO k = nzb+1, nzt+1
1418	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
1419	)*0.286_wp + l_d_cp ql(k,j,i)
1420	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1421
1422	ENDDO
1423
1424	!
1425	!-- Avoid temperature/humidity jumps at the top of the LES domain by
1426	!-- linear interpolation from nzt+2 to nzt+7
1427	DO k = nzt+2, nzt+7
1428	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
1429	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
1430	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
1431	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1432
1433	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
1434	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1435	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
1436	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1437
1438	ENDDO
1439
1440	!-- Linear interpolate to zw grid
1441	DO k = nzb+2, nzt+8
1442	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
1443	rrtm_tlay(0,k-1)) &
1444	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
1445	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1446	ENDDO
1447
1448
1449	!
1450	!-- Calculate liquid water path and cloud fraction for each column.
1451	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
1452	rrtm_cldfr = 0.0_wp
1453	rrtm_reliq = 0.0_wp
1454	rrtm_cliqwp = 0.0_wp
1455	rrtm_icld = 0
1456
1457	DO k = nzb+1, nzt+1
1458	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
1459	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
1460	* 100.0_wp / g
1461
1462	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
1463	rrtm_cldfr(0,k) = 1.0_wp
1464	IF ( rrtm_icld == 0 ) rrtm_icld = 1
1465
1466	!
1467	!-- Calculate cloud droplet effective radius
1468	IF ( cloud_physics ) THEN
1469	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
1470	* rho_surface &
1471	/ ( 4.0_wp * pi * nc_const * rho_l ) &
1472	)**0.33333333333333_wp &
1473	* EXP( LOG( sigma_gc )**2 )
1474
1475	ELSEIF ( cloud_droplets ) THEN
1476	number_of_particles = prt_count(k,j,i)
1477
1478	IF (number_of_particles <= 0) CYCLE
1479	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1480	s_r2 = 0.0_wp
1481	s_r3 = 0.0_wp
1482
1483	DO n = 1, number_of_particles
1484	IF ( particles(n)%particle_mask ) THEN
1485	s_r2 = s_r2 + particles(n)%radius*2 &
1486	particles(n)%weight_factor
1487	s_r3 = s_r3 + particles(n)%radius*3 &
1488	particles(n)%weight_factor
1489	ENDIF
1490	ENDDO
1491
1492	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
1493
1494	ENDIF
1495
1496	!
1497	!-- Limit effective radius
1498	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
1499	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1500	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1501	ENDIF
1502	ENDIF
1503	ENDDO
1504
1505	!
1506	!-- Set surface temperature
1507	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1508
1509	IF ( lw_radiation ) THEN
1510	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
1511	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
1512	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
1513	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
1514	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
1515	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
1516	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
1517	rrtm_reliq , rrtm_lw_tauaer, &
1518	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
1519	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
1520	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
1521
1522	!
1523	!-- Save fluxes
1524	DO k = nzb, nzt+1
1525	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
1526	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1527	ENDDO
1528
1529	!
1530	!-- Save heating rates (convert from K/d to K/h)
1531	DO k = nzb+1, nzt+1
1532	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k) * d_hours_day
1533	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day
1534	ENDDO
1535
1536	!
1537	!-- Save change in LW heating rate
1538	rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1539
1540	ENDIF
1541
1542	IF ( sw_radiation .AND. sun_up ) THEN
1543	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
1544	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
1545	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
1546	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1547	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1548	0.0_wp , day , solar_constant, rrtm_inflgsw,&
1549	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
1550	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
1551	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
1552	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
1553	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
1554	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
1555
1556	!
1557	!-- Save fluxes
1558	DO k = nzb, nzt+1
1559	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
1560	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1561	ENDDO
1562
1563	!
1564	!-- Save heating rates (convert from K/d to K/s)
1565	DO k = nzb+1, nzt+1
1566	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
1567	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
1568	ENDDO
1569
1570	ENDIF
1571
1572	!
1573	!-- Calculate surface net radiation
1574	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
1575	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1576
1577	ENDDO
1578	ENDDO
1579
1580	CALL exchange_horiz( rad_lw_in, nbgp )
1581	CALL exchange_horiz( rad_lw_out, nbgp )
1582	CALL exchange_horiz( rad_lw_hr, nbgp )
1583	CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1584
1585	CALL exchange_horiz( rad_sw_in, nbgp )
1586	CALL exchange_horiz( rad_sw_out, nbgp )
1587	CALL exchange_horiz( rad_sw_hr, nbgp )
1588	CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1589
1590	CALL exchange_horiz_2d( rad_net )
1591	CALL exchange_horiz_2d( rad_lw_out_change_0 )
1592	#endif
1593
1594	END SUBROUTINE radiation_rrtmg
1595
1596
1597	!------------------------------------------------------------------------------!
1598	! Description:
1599	! ------------
1600	!> Calculate the cosine of the zenith angle (variable is called zenith)
1601	!------------------------------------------------------------------------------!
1602	SUBROUTINE calc_zenith
1603
1604	IMPLICIT NONE
1605
1606	REAL(wp) :: declination, & !< solar declination angle
1607	hour_angle !< solar hour angle
1608	!
1609	!-- Calculate current day and time based on the initial values and simulation
1610	!-- time
1611	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
1612	/ 86400.0_wp ), KIND=iwp)
1613	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1614
1615
1616	!
1617	!-- Calculate solar declination and hour angle
1618	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1619	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1620
1621	!
1622	!-- Calculate cosine of solar zenith angle
1623	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
1624	* COS(hour_angle)
1625	zenith(0) = MAX(0.0_wp,zenith(0))
1626
1627	!
1628	!-- Calculate solar directional vector
1629	IF ( sun_direction ) THEN
1630	!-- Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
1631	sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
1632	!-- Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
1633	sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
1634	* COS(declination) * SIN(lat)
1635	ENDIF
1636
1637	!
1638	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
1639	IF ( zenith(0) > 0.0_wp ) THEN
1640	sun_up = .TRUE.
1641	ELSE
1642	sun_up = .FALSE.
1643	END IF
1644
1645	END SUBROUTINE calc_zenith
1646
1647	#if defined ( __rrtmg ) && defined ( __netcdf )
1648	!------------------------------------------------------------------------------!
1649	! Description:
1650	! ------------
1651	!> Calculates surface albedo components based on Briegleb (1992) and
1652	!> Briegleb et al. (1986)
1653	!------------------------------------------------------------------------------!
1654	SUBROUTINE calc_albedo
1655
1656	IMPLICIT NONE
1657
1658	IF ( sun_up ) THEN
1659	!
1660	!-- Ocean
1661	IF ( albedo_type == 1 ) THEN
1662	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
1663	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
1664	* ( zenith(0) - 0.5_wp ) &
1665	* ( zenith(0) - 1.0_wp )
1666	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1667	!
1668	!-- Snow
1669	ELSEIF ( albedo_type == 16 ) THEN
1670	IF ( zenith(0) < 0.5_wp ) THEN
1671	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
1672	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1673	* zenith(0))) - 1.0_wp
1674	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
1675	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1676	* zenith(0))) - 1.0_wp
1677
1678	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1679	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1680	ELSE
1681	rrtm_aldir(0,:,:) = aldif
1682	rrtm_asdir(0,:,:) = asdif
1683	ENDIF
1684	!
1685	!-- Sea ice
1686	ELSEIF ( albedo_type == 15 ) THEN
1687	rrtm_aldir(0,:,:) = aldif
1688	rrtm_asdir(0,:,:) = asdif
1689
1690	!
1691	!-- Asphalt
1692	ELSEIF ( albedo_type == 17 ) THEN
1693	rrtm_aldir(0,:,:) = aldif
1694	rrtm_asdir(0,:,:) = asdif
1695	!
1696	!-- Land surfaces
1697	ELSE
1698	SELECT CASE ( albedo_type )
1699
1700	!
1701	!-- Surface types with strong zenith dependence
1702	CASE ( 1, 2, 3, 4, 11, 12, 13 )
1703	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
1704	(1.0_wp + 0.8_wp * zenith(0))
1705	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
1706	(1.0_wp + 0.8_wp * zenith(0))
1707	!
1708	!-- Surface types with weak zenith dependence
1709	CASE ( 5, 6, 7, 8, 9, 10, 14 )
1710	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
1711	(1.0_wp + 0.2_wp * zenith(0))
1712	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
1713	(1.0_wp + 0.2_wp * zenith(0))
1714
1715	CASE DEFAULT
1716
1717	END SELECT
1718	ENDIF
1719	!
1720	!-- Diffusive albedo is taken from Table 2
1721	rrtm_aldif(0,:,:) = aldif
1722	rrtm_asdif(0,:,:) = asdif
1723
1724	ELSE
1725
1726	rrtm_aldir(0,:,:) = 0.0_wp
1727	rrtm_asdir(0,:,:) = 0.0_wp
1728	rrtm_aldif(0,:,:) = 0.0_wp
1729	rrtm_asdif(0,:,:) = 0.0_wp
1730	ENDIF
1731	END SUBROUTINE calc_albedo
1732
1733	!------------------------------------------------------------------------------!
1734	! Description:
1735	! ------------
1736	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1737	!------------------------------------------------------------------------------!
1738	SUBROUTINE read_sounding_data
1739
1740	IMPLICIT NONE
1741
1742	INTEGER(iwp) :: id, & !< NetCDF id of input file
1743	id_dim_zrad, & !< pressure level id in the NetCDF file
1744	id_var, & !< NetCDF variable id
1745	k, & !< loop index
1746	nz_snd, & !< number of vertical levels in the sounding data
1747	nz_snd_start, & !< start vertical index for sounding data to be used
1748	nz_snd_end !< end vertical index for souding data to be used
1749
1750	REAL(wp) :: t_surface !< actual surface temperature
1751
1752	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1753	t_snd_tmp !< temporary temperature profile (sounding)
1754
1755	!
1756	!-- In case of updates, deallocate arrays first (sufficient to check one
1757	!-- array as the others are automatically allocated). This is required
1758	!-- because nzt_rad might change during the update
1759	IF ( ALLOCATED ( hyp_snd ) ) THEN
1760	DEALLOCATE( hyp_snd )
1761	DEALLOCATE( t_snd )
1762	DEALLOCATE( q_snd )
1763	DEALLOCATE ( rrtm_play )
1764	DEALLOCATE ( rrtm_plev )
1765	DEALLOCATE ( rrtm_tlay )
1766	DEALLOCATE ( rrtm_tlev )
1767
1768	DEALLOCATE ( rrtm_h2ovmr )
1769	DEALLOCATE ( rrtm_cicewp )
1770	DEALLOCATE ( rrtm_cldfr )
1771	DEALLOCATE ( rrtm_cliqwp )
1772	DEALLOCATE ( rrtm_reice )
1773	DEALLOCATE ( rrtm_reliq )
1774	DEALLOCATE ( rrtm_lw_taucld )
1775	DEALLOCATE ( rrtm_lw_tauaer )
1776
1777	DEALLOCATE ( rrtm_lwdflx )
1778	DEALLOCATE ( rrtm_lwdflxc )
1779	DEALLOCATE ( rrtm_lwuflx )
1780	DEALLOCATE ( rrtm_lwuflxc )
1781	DEALLOCATE ( rrtm_lwuflx_dt )
1782	DEALLOCATE ( rrtm_lwuflxc_dt )
1783	DEALLOCATE ( rrtm_lwhr )
1784	DEALLOCATE ( rrtm_lwhrc )
1785
1786	DEALLOCATE ( rrtm_sw_taucld )
1787	DEALLOCATE ( rrtm_sw_ssacld )
1788	DEALLOCATE ( rrtm_sw_asmcld )
1789	DEALLOCATE ( rrtm_sw_fsfcld )
1790	DEALLOCATE ( rrtm_sw_tauaer )
1791	DEALLOCATE ( rrtm_sw_ssaaer )
1792	DEALLOCATE ( rrtm_sw_asmaer )
1793	DEALLOCATE ( rrtm_sw_ecaer )
1794
1795	DEALLOCATE ( rrtm_swdflx )
1796	DEALLOCATE ( rrtm_swdflxc )
1797	DEALLOCATE ( rrtm_swuflx )
1798	DEALLOCATE ( rrtm_swuflxc )
1799	DEALLOCATE ( rrtm_swhr )
1800	DEALLOCATE ( rrtm_swhrc )
1801
1802	ENDIF
1803
1804	!
1805	!-- Open file for reading
1806	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1807	CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1808
1809	!
1810	!-- Inquire dimension of z axis and save in nz_snd
1811	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1812	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1813	CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1814
1815	!
1816	! !-- Allocate temporary array for storing pressure data
1817	ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1818	hyp_snd_tmp = 0.0_wp
1819
1820
1821	!-- Read pressure from file
1822	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1823	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
1824	count = (/nz_snd/) )
1825	CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1826
1827	!
1828	!-- Allocate temporary array for storing temperature data
1829	ALLOCATE( t_snd_tmp(1:nz_snd) )
1830	t_snd_tmp = 0.0_wp
1831
1832	!
1833	!-- Read temperature from file
1834	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1835	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
1836	count = (/nz_snd/) )
1837	CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1838
1839	!
1840	!-- Calculate start of sounding data
1841	nz_snd_start = nz_snd + 1
1842	nz_snd_end = nz_snd + 1
1843
1844	!
1845	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
1846	!-- in Pa, hyp_snd in hPa).
1847	DO k = 1, nz_snd
1848	IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
1849	nz_snd_start = k
1850	EXIT
1851	END IF
1852	END DO
1853
1854	IF ( nz_snd_start <= nz_snd ) THEN
1855	nz_snd_end = nz_snd
1856	END IF
1857
1858
1859	!
1860	!-- Calculate of total grid points for RRTMG calculations
1861	nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1862
1863	!
1864	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
1865	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1866	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1867	ALLOCATE( t_snd(nzb+1:nzt_rad) )
1868	ALLOCATE( q_snd(nzb+1:nzt_rad) )
1869	hyp_snd = 0.0_wp
1870	t_snd = 0.0_wp
1871	q_snd = 0.0_wp
1872
1873	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1874	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1875
1876	nc_stat = NF90_CLOSE( id )
1877
1878	!
1879	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
1880	!-- top of the LES domain. This routine does not consider horizontal or
1881	!-- vertical variability of pressure and temperature
1882	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
1883	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
1884
1885	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1886	DO k = nzb+1, nzt+1
1887	rrtm_play(0,k) = hyp(k) * 0.01_wp
1888	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1889	t_surface )**(1.0_wp/0.286_wp)
1890	ENDDO
1891
1892	DO k = nzt+2, nzt_rad
1893	rrtm_play(0,k) = hyp_snd(k)
1894	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1895	ENDDO
1896	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
1897	1.5 * hyp_snd(nzt_rad) &
1898	- 0.5 * hyp_snd(nzt_rad-1) )
1899	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
1900	0.25_wp * rrtm_plev(0,nzt_rad+1) )
1901
1902	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1903
1904	!
1905	!-- Calculate temperature/humidity levels at top of the LES domain.
1906	!-- Currently, the temperature is taken from sounding data (might lead to a
1907	!-- temperature jump at interface. To do: Humidity is currently not
1908	!-- calculated above the LES domain.
1909	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
1910	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
1911	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1912
1913	DO k = nzt+8, nzt_rad
1914	rrtm_tlay(0,k) = t_snd(k)
1915	rrtm_h2ovmr(0,k) = q_snd(k)
1916	ENDDO
1917	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
1918	- rrtm_tlay(0,nzt_rad-1)
1919	DO k = nzt+9, nzt_rad+1
1920	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
1921	- rrtm_tlay(0,k-1)) &
1922	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
1923	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1924	ENDDO
1925	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1926
1927	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
1928	- rrtm_tlev(0,nzt_rad)
1929	!
1930	!-- Allocate remaining RRTMG arrays
1931	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1932	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1933	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1934	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1935	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1936	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1937	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1938	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1939	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1940	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1941	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1942	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1943	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1944	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1945	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
1946
1947	!
1948	!-- The ice phase is currently not considered in PALM
1949	rrtm_cicewp = 0.0_wp
1950	rrtm_reice = 0.0_wp
1951
1952	!
1953	!-- Set other parameters (move to NAMELIST parameters in the future)
1954	rrtm_lw_tauaer = 0.0_wp
1955	rrtm_lw_taucld = 0.0_wp
1956	rrtm_sw_taucld = 0.0_wp
1957	rrtm_sw_ssacld = 0.0_wp
1958	rrtm_sw_asmcld = 0.0_wp
1959	rrtm_sw_fsfcld = 0.0_wp
1960	rrtm_sw_tauaer = 0.0_wp
1961	rrtm_sw_ssaaer = 0.0_wp
1962	rrtm_sw_asmaer = 0.0_wp
1963	rrtm_sw_ecaer = 0.0_wp
1964
1965
1966	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
1967	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
1968	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
1969	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1970	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1971	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1972
1973	rrtm_swdflx = 0.0_wp
1974	rrtm_swuflx = 0.0_wp
1975	rrtm_swhr = 0.0_wp
1976	rrtm_swuflxc = 0.0_wp
1977	rrtm_swdflxc = 0.0_wp
1978	rrtm_swhrc = 0.0_wp
1979
1980	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
1981	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
1982	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
1983	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1984	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1985	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1986
1987	rrtm_lwdflx = 0.0_wp
1988	rrtm_lwuflx = 0.0_wp
1989	rrtm_lwhr = 0.0_wp
1990	rrtm_lwuflxc = 0.0_wp
1991	rrtm_lwdflxc = 0.0_wp
1992	rrtm_lwhrc = 0.0_wp
1993
1994	ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
1995	ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
1996
1997	rrtm_lwuflx_dt = 0.0_wp
1998	rrtm_lwuflxc_dt = 0.0_wp
1999
2000	END SUBROUTINE read_sounding_data
2001
2002
2003	!------------------------------------------------------------------------------!
2004	! Description:
2005	! ------------
2006	!> Read trace gas data from file
2007	!------------------------------------------------------------------------------!
2008	SUBROUTINE read_trace_gas_data
2009
2010	USE rrsw_ncpar
2011
2012	IMPLICIT NONE
2013
2014	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
2015
2016	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: & !< trace gas names
2017	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
2018	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
2019
2020	INTEGER(iwp) :: id, & !< NetCDF id
2021	k, & !< loop index
2022	m, & !< loop index
2023	n, & !< loop index
2024	nabs, & !< number of absorbers
2025	np, & !< number of pressure levels
2026	id_abs, & !< NetCDF id of the respective absorber
2027	id_dim, & !< NetCDF id of asborber's dimension
2028	id_var !< NetCDf id ot the absorber
2029
2030	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
2031
2032
2033	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
2034	rrtm_play_tmp, & !< temporary array for pressure zu-levels
2035	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
2036	trace_path_tmp !< temporary array for storing trace gas path data
2037
2038	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
2039	trace_mls_path, & !< array for storing trace gas path data
2040	trace_mls_tmp !< temporary array for storing trace gas data
2041
2042
2043	!
2044	!-- In case of updates, deallocate arrays first (sufficient to check one
2045	!-- array as the others are automatically allocated)
2046	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
2047	DEALLOCATE ( rrtm_o3vmr )
2048	DEALLOCATE ( rrtm_co2vmr )
2049	DEALLOCATE ( rrtm_ch4vmr )
2050	DEALLOCATE ( rrtm_n2ovmr )
2051	DEALLOCATE ( rrtm_o2vmr )
2052	DEALLOCATE ( rrtm_cfc11vmr )
2053	DEALLOCATE ( rrtm_cfc12vmr )
2054	DEALLOCATE ( rrtm_cfc22vmr )
2055	DEALLOCATE ( rrtm_ccl4vmr )
2056	ENDIF
2057
2058	!
2059	!-- Allocate trace gas profiles
2060	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
2061	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
2062	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
2063	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
2064	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
2065	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
2066	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
2067	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
2068	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
2069
2070	!
2071	!-- Open file for reading
2072	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
2073	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
2074	!
2075	!-- Inquire dimension ids and dimensions
2076	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
2077	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2078	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
2079	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2080
2081	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
2082	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2083	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
2084	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2085
2086
2087	!
2088	!-- Allocate pressure, and trace gas arrays
2089	ALLOCATE( p_mls(1:np) )
2090	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
2091	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
2092
2093
2094	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
2095	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2096	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
2097	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2098
2099	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
2100	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2101	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
2102	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2103
2104
2105	!
2106	!-- Write absorber amounts (mls) to trace_mls
2107	DO n = 1, num_trace_gases
2108	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
2109
2110	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
2111
2112	!
2113	!-- Replace missing values by zero
2114	WHERE ( trace_mls(n,:) > 2.0_wp )
2115	trace_mls(n,:) = 0.0_wp
2116	END WHERE
2117	END DO
2118
2119	DEALLOCATE ( trace_mls_tmp )
2120
2121	nc_stat = NF90_CLOSE( id )
2122	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
2123
2124	!
2125	!-- Add extra pressure level for calculations of the trace gas paths
2126	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
2127	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
2128
2129	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
2130	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
2131	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
2132	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
2133	* rrtm_plev(0,nzt_rad+1) )
2134
2135	!
2136	!-- Calculate trace gas path (zero at surface) with interpolation to the
2137	!-- sounding levels
2138	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
2139
2140	trace_mls_path(nzb+1,:) = 0.0_wp
2141
2142	DO k = nzb+2, nzt_rad+2
2143	DO m = 1, num_trace_gases
2144	trace_mls_path(k,m) = trace_mls_path(k-1,m)
2145
2146	!
2147	!-- When the pressure level is higher than the trace gas pressure
2148	!-- level, assume that
2149	IF ( rrtm_plev_tmp(k-1) > p_mls(1) ) THEN
2150
2151	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
2152	* ( rrtm_plev_tmp(k-1) &
2153	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
2154	) / g
2155	ENDIF
2156
2157	!
2158	!-- Integrate for each sounding level from the contributing p_mls
2159	!-- levels
2160	DO n = 2, np
2161	!
2162	!-- Limit p_mls so that it is within the model level
2163	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
2164	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
2165	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
2166	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
2167
2168	IF ( p_mls_l > p_mls_u ) THEN
2169
2170	!
2171	!-- Calculate weights for interpolation
2172	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
2173	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
2174	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
2175
2176	!
2177	!-- Add level to trace gas path
2178	trace_mls_path(k,m) = trace_mls_path(k,m) &
2179	+ ( p_wgt_u * trace_mls(m,n) &
2180	+ p_wgt_l * trace_mls(m,n-1) ) &
2181	* (p_mls_l - p_mls_u) / g
2182	ENDIF
2183	ENDDO
2184
2185	IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN
2186	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
2187	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
2188	- rrtm_plev_tmp(k) &
2189	) / g
2190	ENDIF
2191	ENDDO
2192	ENDDO
2193
2194
2195	!
2196	!-- Prepare trace gas path profiles
2197	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
2198
2199	DO m = 1, num_trace_gases
2200
2201	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
2202	- trace_mls_path(1:nzt_rad+1,m) ) * g &
2203	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
2204	- rrtm_plev_tmp(2:nzt_rad+2) )
2205
2206	!
2207	!-- Save trace gas paths to the respective arrays
2208	SELECT CASE ( TRIM( trace_names(m) ) )
2209
2210	CASE ( 'O3' )
2211
2212	rrtm_o3vmr(0,:) = trace_path_tmp(:)
2213
2214	CASE ( 'CO2' )
2215
2216	rrtm_co2vmr(0,:) = trace_path_tmp(:)
2217
2218	CASE ( 'CH4' )
2219
2220	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2221
2222	CASE ( 'N2O' )
2223
2224	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2225
2226	CASE ( 'O2' )
2227
2228	rrtm_o2vmr(0,:) = trace_path_tmp(:)
2229
2230	CASE ( 'CFC11' )
2231
2232	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2233
2234	CASE ( 'CFC12' )
2235
2236	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2237
2238	CASE ( 'CFC22' )
2239
2240	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2241
2242	CASE ( 'CCL4' )
2243
2244	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2245
2246	CASE DEFAULT
2247
2248	END SELECT
2249
2250	ENDDO
2251
2252	DEALLOCATE ( trace_path_tmp )
2253	DEALLOCATE ( trace_mls_path )
2254	DEALLOCATE ( rrtm_play_tmp )
2255	DEALLOCATE ( rrtm_plev_tmp )
2256	DEALLOCATE ( trace_mls )
2257	DEALLOCATE ( p_mls )
2258
2259	END SUBROUTINE read_trace_gas_data
2260
2261
2262	SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2263
2264	USE control_parameters, &
2265	ONLY: message_string
2266
2267	USE NETCDF
2268
2269	USE pegrid
2270
2271	IMPLICIT NONE
2272
2273	CHARACTER(LEN=6) :: message_identifier
2274	CHARACTER(LEN=*) :: routine_name
2275
2276	INTEGER(iwp) :: errno
2277
2278	IF ( nc_stat /= NF90_NOERR ) THEN
2279
2280	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
2281	message_string = TRIM( NF90_STRERROR( nc_stat ) )
2282
2283	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2284
2285	ENDIF
2286
2287	END SUBROUTINE netcdf_handle_error_rad
2288	#endif
2289
2290
2291	!------------------------------------------------------------------------------!
2292	! Description:
2293	! ------------
2294	!> Calculate temperature tendency due to radiative cooling/heating.
2295	!> Cache-optimized version.
2296	!------------------------------------------------------------------------------!
2297	SUBROUTINE radiation_tendency_ij ( i, j, tend )
2298
2299	USE cloud_parameters, &
2300	ONLY: pt_d_t
2301
2302	IMPLICIT NONE
2303
2304	INTEGER(iwp) :: i, j, k !< loop indices
2305
2306	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2307
2308	IF ( radiation_scheme == 'rrtmg' ) THEN
2309	#if defined ( __rrtmg )
2310	!
2311	!-- Calculate tendency based on heating rate
2312	DO k = nzb+1, nzt+1
2313	tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) &
2314	* pt_d_t(k) * d_seconds_hour
2315	ENDDO
2316	#endif
2317	ENDIF
2318
2319	END SUBROUTINE radiation_tendency_ij
2320
2321
2322	!------------------------------------------------------------------------------!
2323	! Description:
2324	! ------------
2325	!> Calculate temperature tendency due to radiative cooling/heating.
2326	!> Vector-optimized version
2327	!------------------------------------------------------------------------------!
2328	SUBROUTINE radiation_tendency ( tend )
2329
2330	USE cloud_parameters, &
2331	ONLY: pt_d_t
2332
2333	USE indices, &
2334	ONLY: nxl, nxr, nyn, nys
2335
2336	IMPLICIT NONE
2337
2338	INTEGER(iwp) :: i, j, k !< loop indices
2339
2340	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2341
2342	IF ( radiation_scheme == 'rrtmg' ) THEN
2343	#if defined ( __rrtmg )
2344	!
2345	!-- Calculate tendency based on heating rate
2346	DO i = nxl, nxr
2347	DO j = nys, nyn
2348	DO k = nzb+1, nzt+1
2349	tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) &
2350	+ rad_sw_hr(k,j,i) ) * pt_d_t(k) &
2351	* d_seconds_hour
2352	ENDDO
2353	ENDDO
2354	ENDDO
2355	#endif
2356	ENDIF
2357
2358
2359	END SUBROUTINE radiation_tendency
2360
2361	!------------------------------------------------------------------------------!
2362	!
2363	! Description:
2364	! ------------
2365	!> Subroutine for averaging 3D data
2366	!------------------------------------------------------------------------------!
2367	SUBROUTINE radiation_3d_data_averaging( mode, variable )
2368
2369
2370	USE control_parameters
2371
2372	USE indices
2373
2374	USE kinds
2375
2376	IMPLICIT NONE
2377
2378	CHARACTER (LEN=*) :: mode !<
2379	CHARACTER (LEN=*) :: variable !<
2380
2381	INTEGER(iwp) :: i !<
2382	INTEGER(iwp) :: j !<
2383	INTEGER(iwp) :: k !<
2384
2385	IF ( mode == 'allocate' ) THEN
2386
2387	SELECT CASE ( TRIM( variable ) )
2388
2389	CASE ( 'rad_net*' )
2390	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
2391	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
2392	ENDIF
2393	rad_net_av = 0.0_wp
2394
2395	CASE ( 'rad_lw_in' )
2396	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
2397	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2398	ENDIF
2399	rad_lw_in_av = 0.0_wp
2400
2401	CASE ( 'rad_lw_out' )
2402	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
2403	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2404	ENDIF
2405	rad_lw_out_av = 0.0_wp
2406
2407	CASE ( 'rad_lw_cs_hr' )
2408	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
2409	ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2410	ENDIF
2411	rad_lw_cs_hr_av = 0.0_wp
2412
2413	CASE ( 'rad_lw_hr' )
2414	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
2415	ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2416	ENDIF
2417	rad_lw_hr_av = 0.0_wp
2418
2419	CASE ( 'rad_sw_in' )
2420	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
2421	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2422	ENDIF
2423	rad_sw_in_av = 0.0_wp
2424
2425	CASE ( 'rad_sw_out' )
2426	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
2427	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2428	ENDIF
2429	rad_sw_out_av = 0.0_wp
2430
2431	CASE ( 'rad_sw_cs_hr' )
2432	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
2433	ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2434	ENDIF
2435	rad_sw_cs_hr_av = 0.0_wp
2436
2437	CASE ( 'rad_sw_hr' )
2438	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
2439	ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2440	ENDIF
2441	rad_sw_hr_av = 0.0_wp
2442
2443	CASE DEFAULT
2444	CONTINUE
2445
2446	END SELECT
2447
2448	ELSEIF ( mode == 'sum' ) THEN
2449
2450	SELECT CASE ( TRIM( variable ) )
2451
2452	CASE ( 'rad_net*' )
2453	DO i = nxlg, nxrg
2454	DO j = nysg, nyng
2455	rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
2456	ENDDO
2457	ENDDO
2458
2459	CASE ( 'rad_lw_in' )
2460	DO i = nxlg, nxrg
2461	DO j = nysg, nyng
2462	DO k = nzb, nzt+1
2463	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
2464	ENDDO
2465	ENDDO
2466	ENDDO
2467
2468	CASE ( 'rad_lw_out' )
2469	DO i = nxlg, nxrg
2470	DO j = nysg, nyng
2471	DO k = nzb, nzt+1
2472	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
2473	ENDDO
2474	ENDDO
2475	ENDDO
2476
2477	CASE ( 'rad_lw_cs_hr' )
2478	DO i = nxlg, nxrg
2479	DO j = nysg, nyng
2480	DO k = nzb, nzt+1
2481	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
2482	ENDDO
2483	ENDDO
2484	ENDDO
2485
2486	CASE ( 'rad_lw_hr' )
2487	DO i = nxlg, nxrg
2488	DO j = nysg, nyng
2489	DO k = nzb, nzt+1
2490	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
2491	ENDDO
2492	ENDDO
2493	ENDDO
2494
2495	CASE ( 'rad_sw_in' )
2496	DO i = nxlg, nxrg
2497	DO j = nysg, nyng
2498	DO k = nzb, nzt+1
2499	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
2500	ENDDO
2501	ENDDO
2502	ENDDO
2503
2504	CASE ( 'rad_sw_out' )
2505	DO i = nxlg, nxrg
2506	DO j = nysg, nyng
2507	DO k = nzb, nzt+1
2508	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
2509	ENDDO
2510	ENDDO
2511	ENDDO
2512
2513	CASE ( 'rad_sw_cs_hr' )
2514	DO i = nxlg, nxrg
2515	DO j = nysg, nyng
2516	DO k = nzb, nzt+1
2517	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
2518	ENDDO
2519	ENDDO
2520	ENDDO
2521
2522	CASE ( 'rad_sw_hr' )
2523	DO i = nxlg, nxrg
2524	DO j = nysg, nyng
2525	DO k = nzb, nzt+1
2526	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
2527	ENDDO
2528	ENDDO
2529	ENDDO
2530
2531	CASE DEFAULT
2532	CONTINUE
2533
2534	END SELECT
2535
2536	ELSEIF ( mode == 'average' ) THEN
2537
2538	SELECT CASE ( TRIM( variable ) )
2539
2540	CASE ( 'rad_net*' )
2541	DO i = nxlg, nxrg
2542	DO j = nysg, nyng
2543	rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
2544	ENDDO
2545	ENDDO
2546
2547	CASE ( 'rad_lw_in' )
2548	DO i = nxlg, nxrg
2549	DO j = nysg, nyng
2550	DO k = nzb, nzt+1
2551	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2552	ENDDO
2553	ENDDO
2554	ENDDO
2555
2556	CASE ( 'rad_lw_out' )
2557	DO i = nxlg, nxrg
2558	DO j = nysg, nyng
2559	DO k = nzb, nzt+1
2560	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2561	ENDDO
2562	ENDDO
2563	ENDDO
2564
2565	CASE ( 'rad_lw_cs_hr' )
2566	DO i = nxlg, nxrg
2567	DO j = nysg, nyng
2568	DO k = nzb, nzt+1
2569	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2570	ENDDO
2571	ENDDO
2572	ENDDO
2573
2574	CASE ( 'rad_lw_hr' )
2575	DO i = nxlg, nxrg
2576	DO j = nysg, nyng
2577	DO k = nzb, nzt+1
2578	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2579	ENDDO
2580	ENDDO
2581	ENDDO
2582
2583	CASE ( 'rad_sw_in' )
2584	DO i = nxlg, nxrg
2585	DO j = nysg, nyng
2586	DO k = nzb, nzt+1
2587	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2588	ENDDO
2589	ENDDO
2590	ENDDO
2591
2592	CASE ( 'rad_sw_out' )
2593	DO i = nxlg, nxrg
2594	DO j = nysg, nyng
2595	DO k = nzb, nzt+1
2596	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2597	ENDDO
2598	ENDDO
2599	ENDDO
2600
2601	CASE ( 'rad_sw_cs_hr' )
2602	DO i = nxlg, nxrg
2603	DO j = nysg, nyng
2604	DO k = nzb, nzt+1
2605	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2606	ENDDO
2607	ENDDO
2608	ENDDO
2609
2610	CASE ( 'rad_sw_hr' )
2611	DO i = nxlg, nxrg
2612	DO j = nysg, nyng
2613	DO k = nzb, nzt+1
2614	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2615	ENDDO
2616	ENDDO
2617	ENDDO
2618
2619	END SELECT
2620
2621	ENDIF
2622
2623	END SUBROUTINE radiation_3d_data_averaging
2624
2625
2626	!------------------------------------------------------------------------------!
2627	!
2628	! Description:
2629	! ------------
2630	!> Subroutine defining appropriate grid for netcdf variables.
2631	!> It is called out from subroutine netcdf.
2632	!------------------------------------------------------------------------------!
2633	SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2634
2635	IMPLICIT NONE
2636
2637	CHARACTER (LEN=*), INTENT(IN) :: var !<
2638	LOGICAL, INTENT(OUT) :: found !<
2639	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
2640	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
2641	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
2642
2643	found = .TRUE.
2644
2645
2646	!
2647	!-- Check for the grid
2648	SELECT CASE ( TRIM( var ) )
2649
2650	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', &
2651	'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', &
2652	'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz', &
2653	'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz', &
2654	'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
2655	grid_x = 'x'
2656	grid_y = 'y'
2657	grid_z = 'zu'
2658
2659	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out', &
2660	'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
2661	'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
2662	'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
2663	grid_x = 'x'
2664	grid_y = 'y'
2665	grid_z = 'zw'
2666
2667
2668	CASE DEFAULT
2669	found = .FALSE.
2670	grid_x = 'none'
2671	grid_y = 'none'
2672	grid_z = 'none'
2673
2674	END SELECT
2675
2676	END SUBROUTINE radiation_define_netcdf_grid
2677
2678	!------------------------------------------------------------------------------!
2679	!
2680	! Description:
2681	! ------------
2682	!> Subroutine defining 3D output variables
2683	!------------------------------------------------------------------------------!
2684	SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, &
2685	local_pf, two_d )
2686
2687	USE indices
2688
2689	USE kinds
2690
2691
2692	IMPLICIT NONE
2693
2694	CHARACTER (LEN=*) :: grid !<
2695	CHARACTER (LEN=*) :: mode !<
2696	CHARACTER (LEN=*) :: variable !<
2697
2698	INTEGER(iwp) :: av !<
2699	INTEGER(iwp) :: i !<
2700	INTEGER(iwp) :: j !<
2701	INTEGER(iwp) :: k !<
2702
2703	LOGICAL :: found !<
2704	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2705
2706	REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
2707
2708	found = .TRUE.
2709
2710	SELECT CASE ( TRIM( variable ) )
2711
2712	CASE ( 'rad_net*_xy' ) ! 2d-array
2713	IF ( av == 0 ) THEN
2714	DO i = nxlg, nxrg
2715	DO j = nysg, nyng
2716	local_pf(i,j,nzb+1) = rad_net(j,i)
2717	ENDDO
2718	ENDDO
2719	ELSE
2720	DO i = nxlg, nxrg
2721	DO j = nysg, nyng
2722	local_pf(i,j,nzb+1) = rad_net_av(j,i)
2723	ENDDO
2724	ENDDO
2725	ENDIF
2726	two_d = .TRUE.
2727	grid = 'zu1'
2728
2729
2730	CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
2731	IF ( av == 0 ) THEN
2732	DO i = nxlg, nxrg
2733	DO j = nysg, nyng
2734	DO k = nzb, nzt+1
2735	local_pf(i,j,k) = rad_lw_in(k,j,i)
2736	ENDDO
2737	ENDDO
2738	ENDDO
2739	ELSE
2740	DO i = nxlg, nxrg
2741	DO j = nysg, nyng
2742	DO k = nzb, nzt+1
2743	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
2744	ENDDO
2745	ENDDO
2746	ENDDO
2747	ENDIF
2748	IF ( mode == 'xy' ) grid = 'zu'
2749
2750	CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
2751	IF ( av == 0 ) THEN
2752	DO i = nxlg, nxrg
2753	DO j = nysg, nyng
2754	DO k = nzb, nzt+1
2755	local_pf(i,j,k) = rad_lw_out(k,j,i)
2756	ENDDO
2757	ENDDO
2758	ENDDO
2759	ELSE
2760	DO i = nxlg, nxrg
2761	DO j = nysg, nyng
2762	DO k = nzb, nzt+1
2763	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
2764	ENDDO
2765	ENDDO
2766	ENDDO
2767	ENDIF
2768	IF ( mode == 'xy' ) grid = 'zu'
2769
2770	CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
2771	IF ( av == 0 ) THEN
2772	DO i = nxlg, nxrg
2773	DO j = nysg, nyng
2774	DO k = nzb, nzt+1
2775	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
2776	ENDDO
2777	ENDDO
2778	ENDDO
2779	ELSE
2780	DO i = nxlg, nxrg
2781	DO j = nysg, nyng
2782	DO k = nzb, nzt+1
2783	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
2784	ENDDO
2785	ENDDO
2786	ENDDO
2787	ENDIF
2788	IF ( mode == 'xy' ) grid = 'zw'
2789
2790	CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
2791	IF ( av == 0 ) THEN
2792	DO i = nxlg, nxrg
2793	DO j = nysg, nyng
2794	DO k = nzb, nzt+1
2795	local_pf(i,j,k) = rad_lw_hr(k,j,i)
2796	ENDDO
2797	ENDDO
2798	ENDDO
2799	ELSE
2800	DO i = nxlg, nxrg
2801	DO j = nysg, nyng
2802	DO k = nzb, nzt+1
2803	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
2804	ENDDO
2805	ENDDO
2806	ENDDO
2807	ENDIF
2808	IF ( mode == 'xy' ) grid = 'zw'
2809
2810	CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
2811	IF ( av == 0 ) THEN
2812	DO i = nxlg, nxrg
2813	DO j = nysg, nyng
2814	DO k = nzb, nzt+1
2815	local_pf(i,j,k) = rad_sw_in(k,j,i)
2816	ENDDO
2817	ENDDO
2818	ENDDO
2819	ELSE
2820	DO i = nxlg, nxrg
2821	DO j = nysg, nyng
2822	DO k = nzb, nzt+1
2823	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2824	ENDDO
2825	ENDDO
2826	ENDDO
2827	ENDIF
2828	IF ( mode == 'xy' ) grid = 'zu'
2829
2830	CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
2831	IF ( av == 0 ) THEN
2832	DO i = nxlg, nxrg
2833	DO j = nysg, nyng
2834	DO k = nzb, nzt+1
2835	local_pf(i,j,k) = rad_sw_out(k,j,i)
2836	ENDDO
2837	ENDDO
2838	ENDDO
2839	ELSE
2840	DO i = nxlg, nxrg
2841	DO j = nysg, nyng
2842	DO k = nzb, nzt+1
2843	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2844	ENDDO
2845	ENDDO
2846	ENDDO
2847	ENDIF
2848	IF ( mode == 'xy' ) grid = 'zu'
2849
2850	CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
2851	IF ( av == 0 ) THEN
2852	DO i = nxlg, nxrg
2853	DO j = nysg, nyng
2854	DO k = nzb, nzt+1
2855	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2856	ENDDO
2857	ENDDO
2858	ENDDO
2859	ELSE
2860	DO i = nxlg, nxrg
2861	DO j = nysg, nyng
2862	DO k = nzb, nzt+1
2863	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2864	ENDDO
2865	ENDDO
2866	ENDDO
2867	ENDIF
2868	IF ( mode == 'xy' ) grid = 'zw'
2869
2870	CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
2871	IF ( av == 0 ) THEN
2872	DO i = nxlg, nxrg
2873	DO j = nysg, nyng
2874	DO k = nzb, nzt+1
2875	local_pf(i,j,k) = rad_sw_hr(k,j,i)
2876	ENDDO
2877	ENDDO
2878	ENDDO
2879	ELSE
2880	DO i = nxlg, nxrg
2881	DO j = nysg, nyng
2882	DO k = nzb, nzt+1
2883	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
2884	ENDDO
2885	ENDDO
2886	ENDDO
2887	ENDIF
2888	IF ( mode == 'xy' ) grid = 'zw'
2889
2890	CASE DEFAULT
2891	found = .FALSE.
2892	grid = 'none'
2893
2894	END SELECT
2895
2896	END SUBROUTINE radiation_data_output_2d
2897
2898
2899	!------------------------------------------------------------------------------!
2900	!
2901	! Description:
2902	! ------------
2903	!> Subroutine defining 3D output variables
2904	!------------------------------------------------------------------------------!
2905	SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
2906
2907
2908	USE indices
2909
2910	USE kinds
2911
2912
2913	IMPLICIT NONE
2914
2915	CHARACTER (LEN=*) :: variable !<
2916
2917	INTEGER(iwp) :: av !<
2918	INTEGER(iwp) :: i !<
2919	INTEGER(iwp) :: j !<
2920	INTEGER(iwp) :: k !<
2921
2922	LOGICAL :: found !<
2923
2924	REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
2925
2926
2927	found = .TRUE.
2928
2929
2930	SELECT CASE ( TRIM( variable ) )
2931
2932	CASE ( 'rad_sw_in' )
2933	IF ( av == 0 ) THEN
2934	DO i = nxlg, nxrg
2935	DO j = nysg, nyng
2936	DO k = nzb, nzt+1
2937	local_pf(i,j,k) = rad_sw_in(k,j,i)
2938	ENDDO
2939	ENDDO
2940	ENDDO
2941	ELSE
2942	DO i = nxlg, nxrg
2943	DO j = nysg, nyng
2944	DO k = nzb, nzt+1
2945	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2946	ENDDO
2947	ENDDO
2948	ENDDO
2949	ENDIF
2950
2951	CASE ( 'rad_sw_out' )
2952	IF ( av == 0 ) THEN
2953	DO i = nxlg, nxrg
2954	DO j = nysg, nyng
2955	DO k = nzb, nzt+1
2956	local_pf(i,j,k) = rad_sw_out(k,j,i)
2957	ENDDO
2958	ENDDO
2959	ENDDO
2960	ELSE
2961	DO i = nxlg, nxrg
2962	DO j = nysg, nyng
2963	DO k = nzb, nzt+1
2964	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2965	ENDDO
2966	ENDDO
2967	ENDDO
2968	ENDIF
2969
2970	CASE ( 'rad_sw_cs_hr' )
2971	IF ( av == 0 ) THEN
2972	DO i = nxlg, nxrg
2973	DO j = nysg, nyng
2974	DO k = nzb, nzt+1
2975	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2976	ENDDO
2977	ENDDO
2978	ENDDO
2979	ELSE
2980	DO i = nxlg, nxrg
2981	DO j = nysg, nyng
2982	DO k = nzb, nzt+1
2983	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2984	ENDDO
2985	ENDDO
2986	ENDDO
2987	ENDIF
2988
2989	CASE ( 'rad_sw_hr' )
2990	IF ( av == 0 ) THEN
2991	DO i = nxlg, nxrg
2992	DO j = nysg, nyng
2993	DO k = nzb, nzt+1
2994	local_pf(i,j,k) = rad_sw_hr(k,j,i)
2995	ENDDO
2996	ENDDO
2997	ENDDO
2998	ELSE
2999	DO i = nxlg, nxrg
3000	DO j = nysg, nyng
3001	DO k = nzb, nzt+1
3002	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3003	ENDDO
3004	ENDDO
3005	ENDDO
3006	ENDIF
3007
3008	CASE ( 'rad_lw_in' )
3009	IF ( av == 0 ) THEN
3010	DO i = nxlg, nxrg
3011	DO j = nysg, nyng
3012	DO k = nzb, nzt+1
3013	local_pf(i,j,k) = rad_lw_in(k,j,i)
3014	ENDDO
3015	ENDDO
3016	ENDDO
3017	ELSE
3018	DO i = nxlg, nxrg
3019	DO j = nysg, nyng
3020	DO k = nzb, nzt+1
3021	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
3022	ENDDO
3023	ENDDO
3024	ENDDO
3025	ENDIF
3026
3027	CASE ( 'rad_lw_out' )
3028	IF ( av == 0 ) THEN
3029	DO i = nxlg, nxrg
3030	DO j = nysg, nyng
3031	DO k = nzb, nzt+1
3032	local_pf(i,j,k) = rad_lw_out(k,j,i)
3033	ENDDO
3034	ENDDO
3035	ENDDO
3036	ELSE
3037	DO i = nxlg, nxrg
3038	DO j = nysg, nyng
3039	DO k = nzb, nzt+1
3040	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
3041	ENDDO
3042	ENDDO
3043	ENDDO
3044	ENDIF
3045
3046	CASE ( 'rad_lw_cs_hr' )
3047	IF ( av == 0 ) THEN
3048	DO i = nxlg, nxrg
3049	DO j = nysg, nyng
3050	DO k = nzb, nzt+1
3051	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
3052	ENDDO
3053	ENDDO
3054	ENDDO
3055	ELSE
3056	DO i = nxlg, nxrg
3057	DO j = nysg, nyng
3058	DO k = nzb, nzt+1
3059	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
3060	ENDDO
3061	ENDDO
3062	ENDDO
3063	ENDIF
3064
3065	CASE ( 'rad_lw_hr' )
3066	IF ( av == 0 ) THEN
3067	DO i = nxlg, nxrg
3068	DO j = nysg, nyng
3069	DO k = nzb, nzt+1
3070	local_pf(i,j,k) = rad_lw_hr(k,j,i)
3071	ENDDO
3072	ENDDO
3073	ENDDO
3074	ELSE
3075	DO i = nxlg, nxrg
3076	DO j = nysg, nyng
3077	DO k = nzb, nzt+1
3078	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
3079	ENDDO
3080	ENDDO
3081	ENDDO
3082	ENDIF
3083
3084	CASE DEFAULT
3085	found = .FALSE.
3086
3087	END SELECT
3088
3089
3090	END SUBROUTINE radiation_data_output_3d
3091
3092	!------------------------------------------------------------------------------!
3093	!
3094	! Description:
3095	! ------------
3096	!> Subroutine defining masked data output
3097	!------------------------------------------------------------------------------!
3098	SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
3099
3100	USE control_parameters
3101
3102	USE indices
3103
3104	USE kinds
3105
3106
3107	IMPLICIT NONE
3108
3109	CHARACTER (LEN=*) :: variable !<
3110
3111	INTEGER(iwp) :: av !<
3112	INTEGER(iwp) :: i !<
3113	INTEGER(iwp) :: j !<
3114	INTEGER(iwp) :: k !<
3115
3116	LOGICAL :: found !<
3117
3118	REAL(wp), &
3119	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
3120	local_pf !<
3121
3122
3123	found = .TRUE.
3124
3125	SELECT CASE ( TRIM( variable ) )
3126
3127
3128	CASE ( 'rad_lw_in' )
3129	IF ( av == 0 ) THEN
3130	DO i = 1, mask_size_l(mid,1)
3131	DO j = 1, mask_size_l(mid,2)
3132	DO k = 1, mask_size_l(mid,3)
3133	local_pf(i,j,k) = rad_lw_in(mask_k(mid,k), &
3134	mask_j(mid,j),mask_i(mid,i))
3135	ENDDO
3136	ENDDO
3137	ENDDO
3138	ELSE
3139	DO i = 1, mask_size_l(mid,1)
3140	DO j = 1, mask_size_l(mid,2)
3141	DO k = 1, mask_size_l(mid,3)
3142	local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k), &
3143	mask_j(mid,j),mask_i(mid,i))
3144	ENDDO
3145	ENDDO
3146	ENDDO
3147	ENDIF
3148
3149	CASE ( 'rad_lw_out' )
3150	IF ( av == 0 ) THEN
3151	DO i = 1, mask_size_l(mid,1)
3152	DO j = 1, mask_size_l(mid,2)
3153	DO k = 1, mask_size_l(mid,3)
3154	local_pf(i,j,k) = rad_lw_out(mask_k(mid,k), &
3155	mask_j(mid,j),mask_i(mid,i))
3156	ENDDO
3157	ENDDO
3158	ENDDO
3159	ELSE
3160	DO i = 1, mask_size_l(mid,1)
3161	DO j = 1, mask_size_l(mid,2)
3162	DO k = 1, mask_size_l(mid,3)
3163	local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k), &
3164	mask_j(mid,j),mask_i(mid,i))
3165	ENDDO
3166	ENDDO
3167	ENDDO
3168	ENDIF
3169
3170	CASE ( 'rad_lw_cs_hr' )
3171	IF ( av == 0 ) THEN
3172	DO i = 1, mask_size_l(mid,1)
3173	DO j = 1, mask_size_l(mid,2)
3174	DO k = 1, mask_size_l(mid,3)
3175	local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k), &
3176	mask_j(mid,j),mask_i(mid,i))
3177	ENDDO
3178	ENDDO
3179	ENDDO
3180	ELSE
3181	DO i = 1, mask_size_l(mid,1)
3182	DO j = 1, mask_size_l(mid,2)
3183	DO k = 1, mask_size_l(mid,3)
3184	local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k), &
3185	mask_j(mid,j),mask_i(mid,i))
3186	ENDDO
3187	ENDDO
3188	ENDDO
3189	ENDIF
3190
3191	CASE ( 'rad_lw_hr' )
3192	IF ( av == 0 ) THEN
3193	DO i = 1, mask_size_l(mid,1)
3194	DO j = 1, mask_size_l(mid,2)
3195	DO k = 1, mask_size_l(mid,3)
3196	local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k), &
3197	mask_j(mid,j),mask_i(mid,i))
3198	ENDDO
3199	ENDDO
3200	ENDDO
3201	ELSE
3202	DO i = 1, mask_size_l(mid,1)
3203	DO j = 1, mask_size_l(mid,2)
3204	DO k = 1, mask_size_l(mid,3)
3205	local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k), &
3206	mask_j(mid,j),mask_i(mid,i))
3207	ENDDO
3208	ENDDO
3209	ENDDO
3210	ENDIF
3211
3212	CASE ( 'rad_sw_in' )
3213	IF ( av == 0 ) THEN
3214	DO i = 1, mask_size_l(mid,1)
3215	DO j = 1, mask_size_l(mid,2)
3216	DO k = 1, mask_size_l(mid,3)
3217	local_pf(i,j,k) = rad_sw_in(mask_k(mid,k), &
3218	mask_j(mid,j),mask_i(mid,i))
3219	ENDDO
3220	ENDDO
3221	ENDDO
3222	ELSE
3223	DO i = 1, mask_size_l(mid,1)
3224	DO j = 1, mask_size_l(mid,2)
3225	DO k = 1, mask_size_l(mid,3)
3226	local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k), &
3227	mask_j(mid,j),mask_i(mid,i))
3228	ENDDO
3229	ENDDO
3230	ENDDO
3231	ENDIF
3232
3233	CASE ( 'rad_sw_out' )
3234	IF ( av == 0 ) THEN
3235	DO i = 1, mask_size_l(mid,1)
3236	DO j = 1, mask_size_l(mid,2)
3237	DO k = 1, mask_size_l(mid,3)
3238	local_pf(i,j,k) = rad_sw_out(mask_k(mid,k), &
3239	mask_j(mid,j),mask_i(mid,i))
3240	ENDDO
3241	ENDDO
3242	ENDDO
3243	ELSE
3244	DO i = 1, mask_size_l(mid,1)
3245	DO j = 1, mask_size_l(mid,2)
3246	DO k = 1, mask_size_l(mid,3)
3247	local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k), &
3248	mask_j(mid,j),mask_i(mid,i))
3249	ENDDO
3250	ENDDO
3251	ENDDO
3252	ENDIF
3253
3254	CASE ( 'rad_sw_cs_hr' )
3255	IF ( av == 0 ) THEN
3256	DO i = 1, mask_size_l(mid,1)
3257	DO j = 1, mask_size_l(mid,2)
3258	DO k = 1, mask_size_l(mid,3)
3259	local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k), &
3260	mask_j(mid,j),mask_i(mid,i))
3261	ENDDO
3262	ENDDO
3263	ENDDO
3264	ELSE
3265	DO i = 1, mask_size_l(mid,1)
3266	DO j = 1, mask_size_l(mid,2)
3267	DO k = 1, mask_size_l(mid,3)
3268	local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k), &
3269	mask_j(mid,j),mask_i(mid,i))
3270	ENDDO
3271	ENDDO
3272	ENDDO
3273	ENDIF
3274
3275	CASE ( 'rad_sw_hr' )
3276	IF ( av == 0 ) THEN
3277	DO i = 1, mask_size_l(mid,1)
3278	DO j = 1, mask_size_l(mid,2)
3279	DO k = 1, mask_size_l(mid,3)
3280	local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k), &
3281	mask_j(mid,j),mask_i(mid,i))
3282	ENDDO
3283	ENDDO
3284	ENDDO
3285	ELSE
3286	DO i = 1, mask_size_l(mid,1)
3287	DO j = 1, mask_size_l(mid,2)
3288	DO k = 1, mask_size_l(mid,3)
3289	local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k), &
3290	mask_j(mid,j),mask_i(mid,i))
3291	ENDDO
3292	ENDDO
3293	ENDDO
3294	ENDIF
3295
3296	CASE DEFAULT
3297	found = .FALSE.
3298
3299	END SELECT
3300
3301
3302	END SUBROUTINE radiation_data_output_mask
3303
3304
3305	!------------------------------------------------------------------------------!
3306	!
3307	! Description:
3308	! ------------
3309	!> Subroutine defines masked output variables
3310	!------------------------------------------------------------------------------!
3311	SUBROUTINE radiation_last_actions
3312
3313
3314	USE control_parameters
3315
3316	USE kinds
3317
3318	IMPLICIT NONE
3319
3320	IF ( write_binary(1:4) == 'true' ) THEN
3321	IF ( ALLOCATED( rad_net ) ) THEN
3322	WRITE ( 14 ) 'rad_net '; WRITE ( 14 ) rad_net
3323	ENDIF
3324	IF ( ALLOCATED( rad_net_av ) ) THEN
3325	WRITE ( 14 ) 'rad_net_av '; WRITE ( 14 ) rad_net_av
3326	ENDIF
3327	IF ( ALLOCATED( rad_lw_in ) ) THEN
3328	WRITE ( 14 ) 'rad_lw_in '; WRITE ( 14 ) rad_lw_in
3329	ENDIF
3330	IF ( ALLOCATED( rad_lw_in_av ) ) THEN
3331	WRITE ( 14 ) 'rad_lw_in_av '; WRITE ( 14 ) rad_lw_in_av
3332	ENDIF
3333	IF ( ALLOCATED( rad_lw_out ) ) THEN
3334	WRITE ( 14 ) 'rad_lw_out '; WRITE ( 14 ) rad_lw_out
3335	ENDIF
3336	IF ( ALLOCATED( rad_lw_out_av ) ) THEN
3337	WRITE ( 14 ) 'rad_lw_out_av '; WRITE ( 14 ) rad_lw_out_av
3338	ENDIF
3339	IF ( ALLOCATED( rad_lw_out_change_0 ) ) THEN
3340	WRITE ( 14 ) 'rad_lw_out_change_0 '
3341	WRITE ( 14 ) rad_lw_out_change_0
3342	ENDIF
3343	IF ( ALLOCATED( rad_lw_cs_hr ) ) THEN
3344	WRITE ( 14 ) 'rad_lw_cs_hr '; WRITE ( 14 ) rad_lw_cs_hr
3345	ENDIF
3346	IF ( ALLOCATED( rad_lw_cs_hr_av ) ) THEN
3347	WRITE ( 14 ) 'rad_lw_cs_hr_av '; WRITE ( 14 ) rad_lw_cs_hr_av
3348	ENDIF
3349	IF ( ALLOCATED( rad_lw_hr ) ) THEN
3350	WRITE ( 14 ) 'rad_lw_hr '; WRITE ( 14 ) rad_lw_hr
3351	ENDIF
3352	IF ( ALLOCATED( rad_lw_hr_av ) ) THEN
3353	WRITE ( 14 ) 'rad_lw_hr_av '; WRITE ( 14 ) rad_lw_hr_av
3354	ENDIF
3355	IF ( ALLOCATED( rad_sw_in ) ) THEN
3356	WRITE ( 14 ) 'rad_sw_in '; WRITE ( 14 ) rad_sw_in
3357	ENDIF
3358	IF ( ALLOCATED( rad_sw_in_av ) ) THEN
3359	WRITE ( 14 ) 'rad_sw_in_av '; WRITE ( 14 ) rad_sw_in_av
3360	ENDIF
3361	IF ( ALLOCATED( rad_sw_out ) ) THEN
3362	WRITE ( 14 ) 'rad_sw_out '; WRITE ( 14 ) rad_sw_out
3363	ENDIF
3364	IF ( ALLOCATED( rad_sw_out_av ) ) THEN
3365	WRITE ( 14 ) 'rad_sw_out_av '; WRITE ( 14 ) rad_sw_out_av
3366	ENDIF
3367	IF ( ALLOCATED( rad_sw_cs_hr ) ) THEN
3368	WRITE ( 14 ) 'rad_sw_cs_hr '; WRITE ( 14 ) rad_sw_cs_hr
3369	ENDIF
3370	IF ( ALLOCATED( rad_sw_cs_hr_av ) ) THEN
3371	WRITE ( 14 ) 'rad_sw_cs_hr_av '; WRITE ( 14 ) rad_sw_cs_hr_av
3372	ENDIF
3373	IF ( ALLOCATED( rad_sw_hr ) ) THEN
3374	WRITE ( 14 ) 'rad_sw_hr '; WRITE ( 14 ) rad_sw_hr
3375	ENDIF
3376	IF ( ALLOCATED( rad_sw_hr_av ) ) THEN
3377	WRITE ( 14 ) 'rad_sw_hr_av '; WRITE ( 14 ) rad_sw_hr_av
3378	ENDIF
3379
3380	WRITE ( 14 ) '* end rad * '
3381
3382	ENDIF
3383
3384	END SUBROUTINE radiation_last_actions
3385
3386
3387	SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
3388	nynfa, nyn_on_file, nysfa, nys_on_file, &
3389	offset_xa, offset_ya, overlap_count, &
3390	tmp_2d, tmp_3d )
3391
3392
3393	USE control_parameters
3394
3395	USE indices
3396
3397	USE kinds
3398
3399	USE pegrid
3400
3401	IMPLICIT NONE
3402
3403	CHARACTER (LEN=20) :: field_char !<
3404
3405	INTEGER(iwp) :: i !<
3406	INTEGER(iwp) :: k !<
3407	INTEGER(iwp) :: nxlc !<
3408	INTEGER(iwp) :: nxlf !<
3409	INTEGER(iwp) :: nxl_on_file !<
3410	INTEGER(iwp) :: nxrc !<
3411	INTEGER(iwp) :: nxrf !<
3412	INTEGER(iwp) :: nxr_on_file !<
3413	INTEGER(iwp) :: nync !<
3414	INTEGER(iwp) :: nynf !<
3415	INTEGER(iwp) :: nyn_on_file !<
3416	INTEGER(iwp) :: nysc !<
3417	INTEGER(iwp) :: nysf !<
3418	INTEGER(iwp) :: nys_on_file !<
3419	INTEGER(iwp) :: overlap_count !<
3420
3421	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
3422	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
3423	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
3424	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
3425	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
3426	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
3427
3428	REAL(wp), &
3429	DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3430	tmp_2d !<
3431
3432	REAL(wp), &
3433	DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3434	tmp_3d !<
3435
3436	REAL(wp), &
3437	DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3438	tmp_3d2 !<
3439
3440
3441
3442	IF ( initializing_actions == 'read_restart_data' ) THEN
3443	READ ( 13 ) field_char
3444
3445	DO WHILE ( TRIM( field_char ) /= '* end rad *' )
3446
3447	DO k = 1, overlap_count
3448
3449	nxlf = nxlfa(i,k)
3450	nxlc = nxlfa(i,k) + offset_xa(i,k)
3451	nxrf = nxrfa(i,k)
3452	nxrc = nxrfa(i,k) + offset_xa(i,k)
3453	nysf = nysfa(i,k)
3454	nysc = nysfa(i,k) + offset_ya(i,k)
3455	nynf = nynfa(i,k)
3456	nync = nynfa(i,k) + offset_ya(i,k)
3457
3458
3459	SELECT CASE ( TRIM( field_char ) )
3460
3461	CASE ( 'rad_net' )
3462	IF ( .NOT. ALLOCATED( rad_net ) ) THEN
3463	ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
3464	ENDIF
3465	IF ( k == 1 ) READ ( 13 ) tmp_2d
3466	rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3467	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3468
3469	CASE ( 'rad_net_av' )
3470	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
3471	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
3472	ENDIF
3473	IF ( k == 1 ) READ ( 13 ) tmp_2d
3474	rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3475	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3476	CASE ( 'rad_lw_in' )
3477	IF ( .NOT. ALLOCATED( rad_lw_in ) ) THEN
3478	IF ( radiation_scheme == 'clear-sky' .OR. &
3479	radiation_scheme == 'constant') THEN
3480	ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
3481	ELSE
3482	ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3483	ENDIF
3484	ENDIF
3485	IF ( k == 1 ) THEN
3486	IF ( radiation_scheme == 'clear-sky' .OR. &
3487	radiation_scheme == 'constant') THEN
3488	READ ( 13 ) tmp_3d2
3489	rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3490	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3491	ELSE
3492	READ ( 13 ) tmp_3d
3493	rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3494	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3495	ENDIF
3496	ENDIF
3497
3498	CASE ( 'rad_lw_in_av' )
3499	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
3500	IF ( radiation_scheme == 'clear-sky' .OR. &
3501	radiation_scheme == 'constant') THEN
3502	ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3503	ELSE
3504	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3505	ENDIF
3506	ENDIF
3507	IF ( k == 1 ) THEN
3508	IF ( radiation_scheme == 'clear-sky' .OR. &
3509	radiation_scheme == 'constant') THEN
3510	READ ( 13 ) tmp_3d2
3511	rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3512	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3513	ELSE
3514	READ ( 13 ) tmp_3d
3515	rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3516	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3517	ENDIF
3518	ENDIF
3519
3520	CASE ( 'rad_lw_out' )
3521	IF ( .NOT. ALLOCATED( rad_lw_out ) ) THEN
3522	IF ( radiation_scheme == 'clear-sky' .OR. &
3523	radiation_scheme == 'constant') THEN
3524	ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
3525	ELSE
3526	ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3527	ENDIF
3528	ENDIF
3529	IF ( k == 1 ) THEN
3530	IF ( radiation_scheme == 'clear-sky' .OR. &
3531	radiation_scheme == 'constant') THEN
3532	READ ( 13 ) tmp_3d2
3533	rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3534	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3535	ELSE
3536	READ ( 13 ) tmp_3d
3537	rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3538	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3539	ENDIF
3540	ENDIF
3541
3542	CASE ( 'rad_lw_out_av' )
3543	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
3544	IF ( radiation_scheme == 'clear-sky' .OR. &
3545	radiation_scheme == 'constant') THEN
3546	ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3547	ELSE
3548	ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3549	ENDIF
3550	ENDIF
3551	IF ( k == 1 ) THEN
3552	IF ( radiation_scheme == 'clear-sky' .OR. &
3553	radiation_scheme == 'constant') THEN
3554	READ ( 13 ) tmp_3d2
3555	rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3556	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3557	ELSE
3558	READ ( 13 ) tmp_3d
3559	rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3560	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3561	ENDIF
3562	ENDIF
3563
3564	CASE ( 'rad_lw_out_change_0' )
3565	IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) ) THEN
3566	ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
3567	ENDIF
3568	IF ( k == 1 ) READ ( 13 ) tmp_2d
3569	rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
3570	= tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3571
3572	CASE ( 'rad_lw_cs_hr' )
3573	IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) ) THEN
3574	ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3575	ENDIF
3576	IF ( k == 1 ) READ ( 13 ) tmp_3d
3577	rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3578	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3579
3580	CASE ( 'rad_lw_cs_hr_av' )
3581	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
3582	ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3583	ENDIF
3584	IF ( k == 1 ) READ ( 13 ) tmp_3d
3585	rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3586	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3587
3588	CASE ( 'rad_lw_hr' )
3589	IF ( .NOT. ALLOCATED( rad_lw_hr ) ) THEN
3590	ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3591	ENDIF
3592	IF ( k == 1 ) READ ( 13 ) tmp_3d
3593	rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3594	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3595
3596	CASE ( 'rad_lw_hr_av' )
3597	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
3598	ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3599	ENDIF
3600	IF ( k == 1 ) READ ( 13 ) tmp_3d
3601	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3602	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3603
3604	CASE ( 'rad_sw_in' )
3605	IF ( .NOT. ALLOCATED( rad_sw_in ) ) THEN
3606	IF ( radiation_scheme == 'clear-sky' .OR. &
3607	radiation_scheme == 'constant') THEN
3608	ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
3609	ELSE
3610	ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3611	ENDIF
3612	ENDIF
3613	IF ( k == 1 ) THEN
3614	IF ( radiation_scheme == 'clear-sky' .OR. &
3615	radiation_scheme == 'constant') THEN
3616	READ ( 13 ) tmp_3d2
3617	rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3618	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3619	ELSE
3620	READ ( 13 ) tmp_3d
3621	rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3622	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3623	ENDIF
3624	ENDIF
3625
3626	CASE ( 'rad_sw_in_av' )
3627	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
3628	IF ( radiation_scheme == 'clear-sky' .OR. &
3629	radiation_scheme == 'constant') THEN
3630	ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3631	ELSE
3632	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3633	ENDIF
3634	ENDIF
3635	IF ( k == 1 ) THEN
3636	IF ( radiation_scheme == 'clear-sky' .OR. &
3637	radiation_scheme == 'constant') THEN
3638	READ ( 13 ) tmp_3d2
3639	rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3640	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3641	ELSE
3642	READ ( 13 ) tmp_3d
3643	rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3644	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3645	ENDIF
3646	ENDIF
3647
3648	CASE ( 'rad_sw_out' )
3649	IF ( .NOT. ALLOCATED( rad_sw_out ) ) THEN
3650	IF ( radiation_scheme == 'clear-sky' .OR. &
3651	radiation_scheme == 'constant') THEN
3652	ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
3653	ELSE
3654	ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3655	ENDIF
3656	ENDIF
3657	IF ( k == 1 ) THEN
3658	IF ( radiation_scheme == 'clear-sky' .OR. &
3659	radiation_scheme == 'constant') THEN
3660	READ ( 13 ) tmp_3d2
3661	rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3662	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3663	ELSE
3664	READ ( 13 ) tmp_3d
3665	rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3666	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3667	ENDIF
3668	ENDIF
3669
3670	CASE ( 'rad_sw_out_av' )
3671	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
3672	IF ( radiation_scheme == 'clear-sky' .OR. &
3673	radiation_scheme == 'constant') THEN
3674	ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3675	ELSE
3676	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3677	ENDIF
3678	ENDIF
3679	IF ( k == 1 ) THEN
3680	IF ( radiation_scheme == 'clear-sky' .OR. &
3681	radiation_scheme == 'constant') THEN
3682	READ ( 13 ) tmp_3d2
3683	rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3684	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3685	ELSE
3686	READ ( 13 ) tmp_3d
3687	rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3688	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3689	ENDIF
3690	ENDIF
3691
3692	CASE ( 'rad_sw_cs_hr' )
3693	IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) ) THEN
3694	ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3695	ENDIF
3696	IF ( k == 1 ) READ ( 13 ) tmp_3d
3697	rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3698	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3699
3700	CASE ( 'rad_sw_cs_hr_av' )
3701	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
3702	ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3703	ENDIF
3704	IF ( k == 1 ) READ ( 13 ) tmp_3d
3705	rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3706	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3707
3708	CASE ( 'rad_sw_hr' )
3709	IF ( .NOT. ALLOCATED( rad_sw_hr ) ) THEN
3710	ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3711	ENDIF
3712	IF ( k == 1 ) READ ( 13 ) tmp_3d
3713	rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3714	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3715
3716	CASE ( 'rad_sw_hr_av' )
3717	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
3718	ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3719	ENDIF
3720	IF ( k == 1 ) READ ( 13 ) tmp_3d
3721	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3722	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3723
3724	CASE DEFAULT
3725	WRITE( message_string, * ) 'unknown variable named "', &
3726	TRIM( field_char ), '" found in', &
3727	'&data from prior run on PE ', myid
3728	CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
3729	0 )
3730
3731	END SELECT
3732
3733	ENDDO
3734
3735	READ ( 13 ) field_char
3736
3737	ENDDO
3738	ENDIF
3739
3740	END SUBROUTINE radiation_read_restart_data
3741
3742
3743	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |