Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |