Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model.f90 @ 1772

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |