Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |