Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |