Home

Context Navigation

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

Last change on this file since 1784 was 1784, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 69.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |