Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |