Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |