Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |