Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |