Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |