Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |