Home

Context Navigation

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

Last change on this file since 1595 was 1591, checked in by maronga, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 63.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |