Home

Context Navigation

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

Last change on this file since 1606 was 1606, checked in by maronga, 9 years ago
bugfix: compilation of radiation_model.f90 without netcdf
Property svn:keywords set to `Id`
File size: 64.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |