Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |