Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |