source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 3885

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/chemistry/SOURCE/radiation_model_mod.f902047-3190,​3218-3297
    /palm/branches/forwind/SOURCE/radiation_model_mod.f901564-1913
    /palm/branches/mosaik_M2/radiation_model_mod.f902360-3471
    /palm/branches/palm4u/SOURCE/radiation_model_mod.f902540-2692
    /palm/branches/radiation/SOURCE/radiation_model_mod.f902081-3493
    /palm/branches/rans/SOURCE/radiation_model_mod.f902078-3128
    /palm/branches/resler/SOURCE/radiation_model_mod.f902023-3605
    /palm/branches/salsa/SOURCE/radiation_model_mod.f902503-3460
    /palm/branches/fricke/SOURCE/radiation_model_mod.f90942-977
    /palm/branches/hoffmann/SOURCE/radiation_model_mod.f90989-1052
    /palm/branches/letzel/masked_output/SOURCE/radiation_model_mod.f90296-409
    /palm/branches/suehring/radiation_model_mod.f90423-666
File size: 503.7 KB
Line 
1!> @file radiation_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2015-2019 Institute of Computer Science of the
18!                     Czech Academy of Sciences, Prague
19! Copyright 2015-2019 Czech Technical University in Prague
20! Copyright 1997-2019 Leibniz Universitaet Hannover
21!------------------------------------------------------------------------------!
22!
23! Current revisions:
24! ------------------
25!
26!
27! Former revisions:
28! -----------------
29! $Id: radiation_model_mod.f90 3885 2019-04-11 11:29:34Z kanani $
30! Changes related to global restructuring of location messages and introduction
31! of additional debug messages
32!
33! 3881 2019-04-10 09:31:22Z suehring
34! Output of albedo and emissivity moved from USM, bugfixes in initialization
35! of albedo
36!
37! 3861 2019-04-04 06:27:41Z maronga
38! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0
39!
40! 3859 2019-04-03 20:30:31Z maronga
41! Added some descriptions
42!
43! 3847 2019-04-01 14:51:44Z suehring
44! Implement check for dt_radiation (must be > 0)
45!
46! 3846 2019-04-01 13:55:30Z suehring
47! unused variable removed
48!
49! 3814 2019-03-26 08:40:31Z pavelkrc
50! Change zenith(0:0) and others to scalar.
51! Code review.
52! Rename exported nzu, nzp and related variables due to name conflict
53!
54! 3771 2019-02-28 12:19:33Z raasch
55! rrtmg preprocessor for directives moved/added, save attribute added to temporary
56! pointers to avoid compiler warnings about outlived pointer targets,
57! statement added to avoid compiler warning about unused variable
58!
59! 3769 2019-02-28 10:16:49Z moh.hefny
60! removed unused variables and subroutine radiation_radflux_gridbox
61!
62! 3767 2019-02-27 08:18:02Z raasch
63! unused variable for file index removed from rrd-subroutines parameter list
64!
65! 3760 2019-02-21 18:47:35Z moh.hefny
66! Bugfix: initialized simulated_time before calculating solar position
67! to enable restart option with reading in SVF from file(s).
68!
69! 3754 2019-02-19 17:02:26Z kanani
70! (resler, pavelkrc)
71! Bugfixes: add further required MRT factors to read/write_svf,
72! fix for aggregating view factors to eliminate local noise in reflected
73! irradiance at mutually close surfaces (corners, presence of trees) in the
74! angular discretization scheme.
75!
76! 3752 2019-02-19 09:37:22Z resler
77! added read/write number of MRT factors to the respective routines
78!
79! 3705 2019-01-29 19:56:39Z suehring
80! Make variables that are sampled in virtual measurement module public
81!
82! 3704 2019-01-29 19:51:41Z suehring
83! Some interface calls moved to module_interface + cleanup
84!
85! 3667 2019-01-10 14:26:24Z schwenkel
86! Modified check for rrtmg input files
87!
88! 3655 2019-01-07 16:51:22Z knoop
89! nopointer option removed
90!
91! 3633 2018-12-17 16:17:57Z schwenkel
92! Include check for rrtmg files
93!
94! 3630 2018-12-17 11:04:17Z knoop
95! - fix initialization of date and time after calling zenith
96! - fix a bug in radiation_solar_pos
97!
98! 3616 2018-12-10 09:44:36Z Salim
99! fix manipulation of time variables in radiation_presimulate_solar_pos
100!
101! 3608 2018-12-07 12:59:57Z suehring $
102! Bugfix radiation output
103!
104! 3607 2018-12-07 11:56:58Z suehring
105! Output of radiation-related quantities migrated to radiation_model_mod.
106!
107! 3589 2018-11-30 15:09:51Z suehring
108! Remove erroneous UTF encoding
109!
110! 3572 2018-11-28 11:40:28Z suehring
111! Add filling the short- and longwave radiation flux arrays (e.g. diffuse,
112! direct, reflected, resedual) for all surfaces. This is required to surface
113! outputs in suface_output_mod. (M. Salim)
114!
115! 3571 2018-11-28 09:24:03Z moh.hefny
116! Add an epsilon value to compare values in if statement to fix possible
117! precsion related errors in raytrace routines.
118!
119! 3524 2018-11-14 13:36:44Z raasch
120! missing cpp-directives added
121!
122! 3495 2018-11-06 15:22:17Z kanani
123! Resort control_parameters ONLY list,
124! From branch radiation@3491 moh.hefny:
125! bugfix in calculating the apparent solar positions by updating
126! the simulated time so that the actual time is correct.
127!
128! 3464 2018-10-30 18:08:55Z kanani
129! From branch resler@3462, pavelkrc:
130! add MRT shaping function for human
131!
132! 3449 2018-10-29 19:36:56Z suehring
133! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague)
134!   - Interaction of plant canopy with LW radiation
135!   - Transpiration from resolved plant canopy dependent on radiation
136!     called from RTM
137!
138!
139! 3435 2018-10-26 18:25:44Z gronemeier
140! - workaround: return unit=illegal in check_data_output for certain variables
141!   when check called from init_masks
142! - Use pointer in masked output to reduce code redundancies
143! - Add terrain-following masked output
144!
145! 3424 2018-10-25 07:29:10Z gronemeier
146! bugfix: add rad_lw_in, rad_lw_out, rad_sw_out to radiation_check_data_output
147!
148! 3378 2018-10-19 12:34:59Z kanani
149! merge from radiation branch (r3362) into trunk
150! (moh.hefny):
151! - removed read/write_svf_on_init and read_dist_max_svf (not used anymore)
152! - bugfix nzut > nzpt in calculating maxboxes
153!
154! 3372 2018-10-18 14:03:19Z raasch
155! bugfix: kind type of 2nd argument of mpi_win_allocate changed, misplaced
156!         __parallel directive
157!
158! 3351 2018-10-15 18:40:42Z suehring
159! Do not overwrite values of spectral and broadband albedo during initialization
160! if they are already initialized in the urban-surface model via ASCII input.
161!
162! 3337 2018-10-12 15:17:09Z kanani
163! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague)
164!   added calculation of the MRT inside the RTM module
165!   MRT fluxes are consequently used in the new biometeorology module
166!   for calculation of biological indices (MRT, PET)
167!   Fixes of v. 2.5 and SVN trunk:
168!    - proper initialization of rad_net_l
169!    - make arrays nsurfs and surfstart TARGET to prevent some MPI problems
170!    - initialization of arrays used in MPI one-sided operation as 1-D arrays
171!      to prevent problems with some MPI/compiler combinations
172!    - fix indexing of target displacement in subroutine request_itarget to
173!      consider nzub
174!    - fix LAD dimmension range in PCB calculation
175!    - check ierr in all MPI calls
176!    - use proper per-gridbox sky and diffuse irradiance
177!    - fix shading for reflected irradiance
178!    - clear away the residuals of "atmospheric surfaces" implementation
179!    - fix rounding bug in raytrace_2d introduced in SVN trunk
180! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague)
181!   can use angular discretization for all SVF
182!   (i.e. reflected and emitted radiation in addition to direct and diffuse),
183!   allowing for much better scaling wih high resoltion and/or complex terrain
184! - Unite array grow factors
185! - Fix slightly shifted terrain height in raytrace_2d
186! - Use more efficient MPI_Win_allocate for reverse gridsurf index
187! - Fix random MPI RMA bugs on Intel compilers
188! - Fix approx. double plant canopy sink values for reflected radiation
189! - Fix mostly missing plant canopy sinks for direct radiation
190! - Fix discretization errors for plant canopy sink in diffuse radiation
191! - Fix rounding errors in raytrace_2d
192!
193! 3274 2018-09-24 15:42:55Z knoop
194! Modularization of all bulk cloud physics code components
195!
196! 3272 2018-09-24 10:16:32Z suehring
197! - split direct and diffusion shortwave radiation using RRTMG rather than using
198!   calc_diffusion_radiation, in case of RRTMG
199! - removed the namelist variable split_diffusion_radiation. Now splitting depends
200!   on the choise of radiation radiation scheme
201! - removed calculating the rdiation flux for surfaces at the radiation scheme
202!   in case of using RTM since it will be calculated anyway in the radiation
203!   interaction routine.
204! - set SW radiation flux for surfaces to zero at night in case of no RTM is used
205! - precalculate the unit vector yxdir of ray direction to avoid the temporarly
206!   array allocation during the subroutine call
207! - fixed a bug in calculating the max number of boxes ray can cross in the domain
208!
209! 3264 2018-09-20 13:54:11Z moh.hefny
210! Bugfix in raytrace_2d calls
211!
212! 3248 2018-09-14 09:42:06Z sward
213! Minor formating changes
214!
215! 3246 2018-09-13 15:14:50Z sward
216! Added error handling for input namelist via parin_fail_message
217!
218! 3241 2018-09-12 15:02:00Z raasch
219! unused variables removed or commented
220!
221! 3233 2018-09-07 13:21:24Z schwenkel
222! Adapted for the use of cloud_droplets
223!
224! 3230 2018-09-05 09:29:05Z schwenkel
225! Bugfix in radiation_constant_surf: changed (10.0 - emissivity_urb) to
226! (1.0 - emissivity_urb)
227!
228! 3226 2018-08-31 12:27:09Z suehring
229! Bugfixes in calculation of sky-view factors and canopy-sink factors.
230!
231! 3186 2018-07-30 17:07:14Z suehring
232! Remove print statement
233!
234! 3180 2018-07-27 11:00:56Z suehring
235! Revise concept for calculation of effective radiative temperature and mapping
236! of radiative heating
237!
238! 3175 2018-07-26 14:07:38Z suehring
239! Bugfix for commit 3172
240!
241! 3173 2018-07-26 12:55:23Z suehring
242! Revise output of surface radiation quantities in case of overhanging
243! structures
244!
245! 3172 2018-07-26 12:06:06Z suehring
246! Bugfixes:
247!  - temporal work-around for calculation of effective radiative surface
248!    temperature
249!  - prevent positive solar radiation during nighttime
250!
251! 3170 2018-07-25 15:19:37Z suehring
252! Bugfix, map signle-column radiation forcing profiles on top of any topography
253!
254! 3156 2018-07-19 16:30:54Z knoop
255! Bugfix: replaced usage of the pt array with the surf%pt_surface array
256!
257! 3137 2018-07-17 06:44:21Z maronga
258! String length for trace_names fixed
259!
260! 3127 2018-07-15 08:01:25Z maronga
261! A few pavement parameters updated.
262!
263! 3123 2018-07-12 16:21:53Z suehring
264! Correct working precision for INTEGER number
265!
266! 3122 2018-07-11 21:46:41Z maronga
267! Bugfix: maximum distance for raytracing was set to  -999 m by default,
268! effectively switching off all surface reflections when max_raytracing_dist
269! was not explicitly set in namelist
270!
271! 3117 2018-07-11 09:59:11Z maronga
272! Bugfix: water vapor was not transfered to RRTMG when bulk_cloud_model = .F.
273! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
274! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
275!
276! 3116 2018-07-10 14:31:58Z suehring
277! Output of long/shortwave radiation at surface
278!
279! 3107 2018-07-06 15:55:51Z suehring
280! Bugfix, missing index for dz
281!
282! 3066 2018-06-12 08:55:55Z Giersch
283! Error message revised
284!
285! 3065 2018-06-12 07:03:02Z Giersch
286! dz was replaced by dz(1), error message concerning vertical stretching was
287! added 
288!
289! 3049 2018-05-29 13:52:36Z Giersch
290! Error messages revised
291!
292! 3045 2018-05-28 07:55:41Z Giersch
293! Error message revised
294!
295! 3026 2018-05-22 10:30:53Z schwenkel
296! Changed the name specific humidity to mixing ratio, since we are computing
297! mixing ratios.
298!
299! 3016 2018-05-09 10:53:37Z Giersch
300! Revised structure of reading svf data according to PALM coding standard:
301! svf_code_field/len and fsvf removed, error messages PA0493 and PA0494 added,
302! allocation status of output arrays checked.
303!
304! 3014 2018-05-09 08:42:38Z maronga
305! Introduced plant canopy height similar to urban canopy height to limit
306! the memory requirement to allocate lad.
307! Deactivated automatic setting of minimum raytracing distance.
308!
309! 3004 2018-04-27 12:33:25Z Giersch
310! Further allocation checks implemented (averaged data will be assigned to fill
311! values if no allocation happened so far)
312!
313! 2995 2018-04-19 12:13:16Z Giersch
314! IF-statement in radiation_init removed so that the calculation of radiative
315! fluxes at model start is done in any case, bugfix in
316! radiation_presimulate_solar_pos (end_time is the sum of end_time and the
317! spinup_time specified in the p3d_file ), list of variables/fields that have
318! to be written out or read in case of restarts has been extended
319!
320! 2977 2018-04-17 10:27:57Z kanani
321! Implement changes from branch radiation (r2948-2971) with minor modifications,
322! plus some formatting.
323! (moh.hefny):
324! - replaced plant_canopy by npcbl to check tree existence to avoid weird
325!   allocation of related arrays (after domain decomposition some domains
326!   contains no trees although plant_canopy (global parameter) is still TRUE).
327! - added a namelist parameter to force RTM settings
328! - enabled the option to switch radiation reflections off
329! - renamed surf_reflections to surface_reflections
330! - removed average_radiation flag from the namelist (now it is implicitly set
331!   in init_3d_model according to RTM)
332! - edited read and write sky view factors and CSF routines to account for
333!   the sub-domains which may not contain any of them
334!
335! 2967 2018-04-13 11:22:08Z raasch
336! bugfix: missing parallel cpp-directives added
337!
338! 2964 2018-04-12 16:04:03Z Giersch
339! Error message PA0491 has been introduced which could be previously found in
340! check_open. The variable numprocs_previous_run is only known in case of
341! initializing_actions == read_restart_data
342!
343! 2963 2018-04-12 14:47:44Z suehring
344! - Introduce index for vegetation/wall, pavement/green-wall and water/window
345!   surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
346! - Minor bugfix in initialization of albedo for window surfaces
347!
348! 2944 2018-04-03 16:20:18Z suehring
349! Fixed bad commit
350!
351! 2943 2018-04-03 16:17:10Z suehring
352! No read of nsurfl from SVF file since it is calculated in
353! radiation_interaction_init,
354! allocation of arrays in radiation_read_svf only if not yet allocated,
355! update of 2920 revision comment.
356!
357! 2932 2018-03-26 09:39:22Z maronga
358! renamed radiation_par to radiation_parameters
359!
360! 2930 2018-03-23 16:30:46Z suehring
361! Remove default surfaces from radiation model, does not make much sense to
362! apply radiation model without energy-balance solvers; Further, add check for
363! this.
364!
365! 2920 2018-03-22 11:22:01Z kanani
366! - Bugfix: Initialize pcbl array (=-1)
367! RTM version 2.0 (Jaroslav Resler, Pavel Krc, Mohamed Salim):
368! - new major version of radiation interactions
369! - substantially enhanced performance and scalability
370! - processing of direct and diffuse solar radiation separated from reflected
371!   radiation, removed virtual surfaces
372! - new type of sky discretization by azimuth and elevation angles
373! - diffuse radiation processed cumulatively using sky view factor
374! - used precalculated apparent solar positions for direct irradiance
375! - added new 2D raytracing process for processing whole vertical column at once
376!   to increase memory efficiency and decrease number of MPI RMA operations
377! - enabled limiting the number of view factors between surfaces by the distance
378!   and value
379! - fixing issues induced by transferring radiation interactions from
380!   urban_surface_mod to radiation_mod
381! - bugfixes and other minor enhancements
382!
383! 2906 2018-03-19 08:56:40Z Giersch
384! NAMELIST paramter read/write_svf_on_init have been removed, functions
385! check_open and close_file are used now for opening/closing files related to
386! svf data, adjusted unit number and error numbers
387!
388! 2894 2018-03-15 09:17:58Z Giersch
389! Calculations of the index range of the subdomain on file which overlaps with
390! the current subdomain are already done in read_restart_data_mod
391! radiation_read_restart_data was renamed to radiation_rrd_local and
392! radiation_last_actions was renamed to radiation_wrd_local, variable named
393! found has been introduced for checking if restart data was found, reading
394! of restart strings has been moved completely to read_restart_data_mod,
395! radiation_rrd_local is already inside the overlap loop programmed in
396! read_restart_data_mod, the marker *** end rad *** is not necessary anymore,
397! strings and their respective lengths are written out and read now in case of
398! restart runs to get rid of prescribed character lengths (Giersch)
399!
400! 2809 2018-02-15 09:55:58Z suehring
401! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
402!
403! 2753 2018-01-16 14:16:49Z suehring
404! Tile approach for spectral albedo implemented.
405!
406! 2746 2018-01-15 12:06:04Z suehring
407! Move flag plant canopy to modules
408!
409! 2724 2018-01-05 12:12:38Z maronga
410! Set default of average_radiation to .FALSE.
411!
412! 2723 2018-01-05 09:27:03Z maronga
413! Bugfix in calculation of rad_lw_out (clear-sky). First grid level was used
414! instead of the surface value
415!
416! 2718 2018-01-02 08:49:38Z maronga
417! Corrected "Former revisions" section
418!
419! 2707 2017-12-18 18:34:46Z suehring
420! Changes from last commit documented
421!
422! 2706 2017-12-18 18:33:49Z suehring
423! Bugfix, in average radiation case calculate exner function before using it.
424!
425! 2701 2017-12-15 15:40:50Z suehring
426! Changes from last commit documented
427!
428! 2698 2017-12-14 18:46:24Z suehring
429! Bugfix in get_topography_top_index
430!
431! 2696 2017-12-14 17:12:51Z kanani
432! - Change in file header (GPL part)
433! - Improved reading/writing of SVF from/to file (BM)
434! - Bugfixes concerning RRTMG as well as average_radiation options (M. Salim)
435! - Revised initialization of surface albedo and some minor bugfixes (MS)
436! - Update net radiation after running radiation interaction routine (MS)
437! - Revisions from M Salim included
438! - Adjustment to topography and surface structure (MS)
439! - Initialization of albedo and surface emissivity via input file (MS)
440! - albedo_pars extended (MS)
441!
442! 2604 2017-11-06 13:29:00Z schwenkel
443! bugfix for calculation of effective radius using morrison microphysics
444!
445! 2601 2017-11-02 16:22:46Z scharf
446! added emissivity to namelist
447!
448! 2575 2017-10-24 09:57:58Z maronga
449! Bugfix: calculation of shortwave and longwave albedos for RRTMG swapped
450!
451! 2547 2017-10-16 12:41:56Z schwenkel
452! extended by cloud_droplets option, minor bugfix and correct calculation of
453! cloud droplet number concentration
454!
455! 2544 2017-10-13 18:09:32Z maronga
456! Moved date and time quantitis to separate module date_and_time_mod
457!
458! 2512 2017-10-04 08:26:59Z raasch
459! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
460! no output of ghost layer data
461!
462! 2504 2017-09-27 10:36:13Z maronga
463! Updates pavement types and albedo parameters
464!
465! 2328 2017-08-03 12:34:22Z maronga
466! Emissivity can now be set individually for each pixel.
467! Albedo type can be inferred from land surface model.
468! Added default albedo type for bare soil
469!
470! 2318 2017-07-20 17:27:44Z suehring
471! Get topography top index via Function call
472!
473! 2317 2017-07-20 17:27:19Z suehring
474! Improved syntax layout
475!
476! 2298 2017-06-29 09:28:18Z raasch
477! type of write_binary changed from CHARACTER to LOGICAL
478!
479! 2296 2017-06-28 07:53:56Z maronga
480! Added output of rad_sw_out for radiation_scheme = 'constant'
481!
482! 2270 2017-06-09 12:18:47Z maronga
483! Numbering changed (2 timeseries removed)
484!
485! 2249 2017-06-06 13:58:01Z sward
486! Allow for RRTMG runs without humidity/cloud physics
487!
488! 2248 2017-06-06 13:52:54Z sward
489! Error no changed
490!
491! 2233 2017-05-30 18:08:54Z suehring
492!
493! 2232 2017-05-30 17:47:52Z suehring
494! Adjustments to new topography concept
495! Bugfix in read restart
496!
497! 2200 2017-04-11 11:37:51Z suehring
498! Bugfix in call of exchange_horiz_2d and read restart data
499!
500! 2163 2017-03-01 13:23:15Z schwenkel
501! Bugfix in radiation_check_data_output
502!
503! 2157 2017-02-22 15:10:35Z suehring
504! Bugfix in read_restart data
505!
506! 2011 2016-09-19 17:29:57Z kanani
507! Removed CALL of auxiliary SUBROUTINE get_usm_info,
508! flag urban_surface is now defined in module control_parameters.
509!
510! 2007 2016-08-24 15:47:17Z kanani
511! Added calculation of solar directional vector for new urban surface
512! model,
513! accounted for urban_surface model in radiation_check_parameters,
514! correction of comments for zenith angle.
515!
516! 2000 2016-08-20 18:09:15Z knoop
517! Forced header and separation lines into 80 columns
518!
519! 1976 2016-07-27 13:28:04Z maronga
520! Output of 2D/3D/masked data is now directly done within this module. The
521! radiation schemes have been simplified for better usability so that
522! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
523! the radiation code used.
524!
525! 1856 2016-04-13 12:56:17Z maronga
526! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
527!
528! 1853 2016-04-11 09:00:35Z maronga
529! Added routine for radiation_scheme = constant.
530
531! 1849 2016-04-08 11:33:18Z hoffmann
532! Adapted for modularization of microphysics
533!
534! 1826 2016-04-07 12:01:39Z maronga
535! Further modularization.
536!
537! 1788 2016-03-10 11:01:04Z maronga
538! Added new albedo class for pavements / roads.
539!
540! 1783 2016-03-06 18:36:17Z raasch
541! palm-netcdf-module removed in order to avoid a circular module dependency,
542! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
543! added
544!
545! 1757 2016-02-22 15:49:32Z maronga
546! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
547! profiles for pressure and temperature above the LES domain.
548!
549! 1709 2015-11-04 14:47:01Z maronga
550! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
551! corrections
552!
553! 1701 2015-11-02 07:43:04Z maronga
554! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
555!
556! 1691 2015-10-26 16:17:44Z maronga
557! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
558! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
559! Added output of radiative heating rates.
560!
561! 1682 2015-10-07 23:56:08Z knoop
562! Code annotations made doxygen readable
563!
564! 1606 2015-06-29 10:43:37Z maronga
565! Added preprocessor directive __netcdf to allow for compiling without netCDF.
566! Note, however, that RRTMG cannot be used without netCDF.
567!
568! 1590 2015-05-08 13:56:27Z maronga
569! Bugfix: definition of character strings requires same length for all elements
570!
571! 1587 2015-05-04 14:19:01Z maronga
572! Added albedo class for snow
573!
574! 1585 2015-04-30 07:05:52Z maronga
575! Added support for RRTMG
576!
577! 1571 2015-03-12 16:12:49Z maronga
578! Added missing KIND attribute. Removed upper-case variable names
579!
580! 1551 2015-03-03 14:18:16Z maronga
581! Added support for data output. Various variables have been renamed. Added
582! interface for different radiation schemes (currently: clear-sky, constant, and
583! RRTM (not yet implemented).
584!
585! 1496 2014-12-02 17:25:50Z maronga
586! Initial revision
587!
588!
589! Description:
590! ------------
591!> Radiation models and interfaces
592!> @todo Replace dz(1) appropriatly to account for grid stretching
593!> @todo move variable definitions used in radiation_init only to the subroutine
594!>       as they are no longer required after initialization.
595!> @todo Output of full column vertical profiles used in RRTMG
596!> @todo Output of other rrtm arrays (such as volume mixing ratios)
597!> @todo Check for mis-used NINT() calls in raytrace_2d
598!>       RESULT: Original was correct (carefully verified formula), the change
599!>               to INT broke raytracing      -- P. Krc
600!> @todo Optimize radiation_tendency routines
601!>
602!> @note Many variables have a leading dummy dimension (0:0) in order to
603!>       match the assume-size shape expected by the RRTMG model.
604!------------------------------------------------------------------------------!
605 MODULE radiation_model_mod
606 
607    USE arrays_3d,                                                             &
608        ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
609
610    USE basic_constants_and_equations_mod,                                     &
611        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
612               barometric_formula
613
614    USE calc_mean_profile_mod,                                                 &
615        ONLY:  calc_mean_profile
616
617    USE control_parameters,                                                    &
618        ONLY:  cloud_droplets, coupling_char,                                  &
619               debug_output, debug_string,                                     &
620               dz, dt_spinup, end_time,                                        &
621               humidity,                                                       &
622               initializing_actions, io_blocks, io_group,                      &
623               land_surface, large_scale_forcing,                              &
624               latitude, longitude, lsf_surf,                                  &
625               message_string, plant_canopy, pt_surface,                       &
626               rho_surface, simulated_time, spinup_time, surface_pressure,     &
627               read_svf, write_svf,                                            &
628               time_since_reference_point, urban_surface, varnamelength
629
630    USE cpulog,                                                                &
631        ONLY:  cpu_log, log_point, log_point_s
632
633    USE grid_variables,                                                        &
634         ONLY:  ddx, ddy, dx, dy
635
636    USE date_and_time_mod,                                                     &
637        ONLY:  calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,&
638               day_of_year, d_seconds_year, day_of_month, day_of_year_init,    &
639               init_date_and_time, month_of_year, time_utc_init, time_utc
640
641    USE indices,                                                               &
642        ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
643               nzb, nzt
644
645    USE, INTRINSIC :: iso_c_binding
646
647    USE kinds
648
649    USE bulk_cloud_model_mod,                                                  &
650        ONLY:  bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
651
652#if defined ( __netcdf )
653    USE NETCDF
654#endif
655
656    USE netcdf_data_input_mod,                                                 &
657        ONLY:  albedo_type_f, albedo_pars_f, building_type_f, pavement_type_f, &
658               vegetation_type_f, water_type_f
659
660    USE plant_canopy_model_mod,                                                &
661        ONLY:  lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate,  &
662               plant_canopy_transpiration, pcm_calc_transpiration_rate
663
664    USE pegrid
665
666#if defined ( __rrtmg )
667    USE parrrsw,                                                               &
668        ONLY:  naerec, nbndsw
669
670    USE parrrtm,                                                               &
671        ONLY:  nbndlw
672
673    USE rrtmg_lw_init,                                                         &
674        ONLY:  rrtmg_lw_ini
675
676    USE rrtmg_sw_init,                                                         &
677        ONLY:  rrtmg_sw_ini
678
679    USE rrtmg_lw_rad,                                                          &
680        ONLY:  rrtmg_lw
681
682    USE rrtmg_sw_rad,                                                          &
683        ONLY:  rrtmg_sw
684#endif
685    USE statistics,                                                            &
686        ONLY:  hom
687
688    USE surface_mod,                                                           &
689        ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
690               ind_pav_green, ind_veg_wall, ind_wat_win,                       &
691               surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
692               vertical_surfaces_exist
693
694    IMPLICIT NONE
695
696    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
697
698!
699!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
700    CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
701                                   'user defined                         ', & !  0
702                                   'ocean                                ', & !  1
703                                   'mixed farming, tall grassland        ', & !  2
704                                   'tall/medium grassland                ', & !  3
705                                   'evergreen shrubland                  ', & !  4
706                                   'short grassland/meadow/shrubland     ', & !  5
707                                   'evergreen needleleaf forest          ', & !  6
708                                   'mixed deciduous evergreen forest     ', & !  7
709                                   'deciduous forest                     ', & !  8
710                                   'tropical evergreen broadleaved forest', & !  9
711                                   'medium/tall grassland/woodland       ', & ! 10
712                                   'desert, sandy                        ', & ! 11
713                                   'desert, rocky                        ', & ! 12
714                                   'tundra                               ', & ! 13
715                                   'land ice                             ', & ! 14
716                                   'sea ice                              ', & ! 15
717                                   'snow                                 ', & ! 16
718                                   'bare soil                            ', & ! 17
719                                   'asphalt/concrete mix                 ', & ! 18
720                                   'asphalt (asphalt concrete)           ', & ! 19
721                                   'concrete (Portland concrete)         ', & ! 20
722                                   'sett                                 ', & ! 21
723                                   'paving stones                        ', & ! 22
724                                   'cobblestone                          ', & ! 23
725                                   'metal                                ', & ! 24
726                                   'wood                                 ', & ! 25
727                                   'gravel                               ', & ! 26
728                                   'fine gravel                          ', & ! 27
729                                   'pebblestone                          ', & ! 28
730                                   'woodchips                            ', & ! 29
731                                   'tartan (sports)                      ', & ! 30
732                                   'artifical turf (sports)              ', & ! 31
733                                   'clay (sports)                        ', & ! 32
734                                   'building (dummy)                     '  & ! 33
735                                                         /)
736
737    INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
738                    dots_rad     = 0_iwp              !< starting index for timeseries output
739
740    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
741                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
742                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
743                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
744                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
745                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
746                sw_radiation = .TRUE.,                & !< flag parameter indicating whether shortwave radiation shall be calculated
747                sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
748                average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
749                radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
750                surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
751                                                        !< When it switched off, only the effect of buildings and trees shadow
752                                                        !< will be considered. However fewer SVFs are expected.
753                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
754
755    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
756                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
757                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
758                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
759                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
760                decl_1,                          & !< declination coef. 1
761                decl_2,                          & !< declination coef. 2
762                decl_3,                          & !< declination coef. 3
763                dt_radiation = 0.0_wp,           & !< radiation model timestep
764                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
765                lon = 0.0_wp,                    & !< longitude in radians
766                lat = 0.0_wp,                    & !< latitude in radians
767                net_radiation = 0.0_wp,          & !< net radiation at surface
768                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
769                sky_trans,                       & !< sky transmissivity
770                time_radiation = 0.0_wp            !< time since last call of radiation code
771
772
773    REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
774    REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
775    REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector
776
777    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
778    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av  !< average of incoming longwave radiation at surface
779    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
780    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
781    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
782
783    REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
784!
785!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)     
786!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
787    REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 
788                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
789                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
790                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
791                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
792                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
793                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
794                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
795                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
796                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
797                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
798                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
799                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
800                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
801                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
802                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
803                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
804                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
805                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 18
806                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 19
807                                   0.30_wp, 0.30_wp, 0.30_wp,            & ! 20
808                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 21
809                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 22
810                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 23
811                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 24
812                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 25
813                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 26
814                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 27
815                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 28
816                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 29
817                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 30
818                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 31
819                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 32
820                                   0.17_wp, 0.17_wp, 0.17_wp             & ! 33
821                                 /), (/ 3, 33 /) )
822
823    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
824                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
825                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
826                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
827                        rad_lw_hr_av,                  & !< average of rad_sw_hr
828                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
829                        rad_lw_in_av,                  & !< average of rad_lw_in
830                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
831                        rad_lw_out_av,                 & !< average of rad_lw_out
832                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
833                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
834                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
835                        rad_sw_hr_av,                  & !< average of rad_sw_hr
836                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
837                        rad_sw_in_av,                  & !< average of rad_sw_in
838                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
839                        rad_sw_out_av                    !< average of rad_sw_out
840
841
842!
843!-- Variables and parameters used in RRTMG only
844#if defined ( __rrtmg )
845    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
846
847
848!
849!-- Flag parameters to be passed to RRTMG (should not be changed until ice phase in clouds is allowed)
850    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
851                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
852                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
853                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
854                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
855                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
856                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
857
858!
859!-- The following variables should be only changed with care, as this will
860!-- require further setting of some variables, which is currently not
861!-- implemented (aerosols, ice phase).
862    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
863                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
864                    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)
865
866    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
867
868    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
869    LOGICAL :: sw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg sw file exists
870    LOGICAL :: lw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg lw file exists
871
872
873    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
874
875    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
876                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
877                                           t_snd          !< actual temperature from sounding data (hPa)
878
879    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
880                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
881                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
882                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
883                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
884                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m2)
885                                             rrtm_cldfr,     & !< cloud fraction (0,1)
886                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m2)
887                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
888                                             rrtm_emis,      & !< surface emissivity (0-1) 
889                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
890                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
891                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
892                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
893                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
894                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
895                                             rrtm_reice,     & !< cloud ice effective radius (microns)
896                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
897                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
898                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
899                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
900                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
901                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
902                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
903                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
904                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
905                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
906                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
907                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
908                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
909                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
910                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
911                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
912                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
913                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2)
914                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2)
915
916    REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
917                                             rrtm_aldir,     & !< surface albedo for longwave direct radiation
918                                             rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
919                                             rrtm_asdir        !< surface albedo for shortwave direct radiation
920
921!
922!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
923    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
924                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
925                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
926                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
927                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
928                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
929                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
930                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
931                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
932                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
933                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
934                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
935                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
936                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
937
938#endif
939!
940!-- Parameters of urban and land surface models
941    INTEGER(iwp)                                   ::  nz_urban                           !< number of layers of urban surface (will be calculated)
942    INTEGER(iwp)                                   ::  nz_plant                           !< number of layers of plant canopy (will be calculated)
943    INTEGER(iwp)                                   ::  nz_urban_b                         !< bottom layer of urban surface (will be calculated)
944    INTEGER(iwp)                                   ::  nz_urban_t                         !< top layer of urban surface (will be calculated)
945    INTEGER(iwp)                                   ::  nz_plant_t                         !< top layer of plant canopy (will be calculated)
946!-- parameters of urban and land surface models
947    INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
948    INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
949    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
950    INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
951    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
952    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
953    INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
954    INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
955    INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
956    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
957    INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
958    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
959
960    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
961
962    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
963    INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
964    INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
965    INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
966    INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
967    INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
968
969    INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
970    INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
971    INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
972    INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
973    INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
974
975    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
976    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
977    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
978    REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
979                                                                                          !< direction (will be calc'd)
980
981
982!-- indices and sizes of urban and land surface models
983    INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
984    INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
985    INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
986    INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
987    INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
988    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
989
990!-- indices needed for RTM netcdf output subroutines
991    INTEGER(iwp), PARAMETER                        :: nd = 5
992    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
993    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
994    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
995    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
996    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
997
998!-- indices and sizes of urban and land surface models
999    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
1000    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear     !< dtto (linearly allocated array)
1001    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
1002    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear      !< dtto (linearly allocated array)
1003    INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
1004    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
1005    INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
1006    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
1007                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
1008
1009!-- block variables needed for calculation of the plant canopy model inside the urban surface model
1010    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
1011    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
1012    INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
1013    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
1014    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
1015    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir       !< array of absorbed direct sw radiation for local plant canopy box
1016    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif       !< array of absorbed diffusion sw radiation for local plant canopy box
1017    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
1018
1019!-- configuration parameters (they can be setup in PALM config)
1020    LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
1021    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
1022                                                                                          !< reflected radiation (as opposed to all mutually visible pairs)
1023    LOGICAL                                        ::  plant_lw_interact = .TRUE.         !< whether plant canopy interacts with LW radiation (in addition to SW)
1024    INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
1025    LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
1026    LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
1027    LOGICAL                                        ::  mrt_geom_human = .TRUE.            !< MRT direction weights simulate human body instead of a sphere
1028    INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
1029    REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
1030    INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
1031    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
1032    INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
1033    INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
1034    REAL(wp)                                       ::  max_raytracing_dist = -999.0_wp    !< maximum distance for raytracing (in metres)
1035    REAL(wp)                                       ::  min_irrf_value = 1e-6_wp           !< minimum potential irradiance factor value for raytracing
1036    REAL(wp), DIMENSION(1:30)                      ::  svfnorm_report_thresh = 1e21_wp    !< thresholds of SVF normalization values to report
1037    INTEGER(iwp)                                   ::  svfnorm_report_num                 !< number of SVF normalization thresholds to report
1038
1039!-- radiation related arrays to be used in radiation_interaction routine
1040    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_dir    !< direct sw radiation
1041    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_diff   !< diffusion sw radiation
1042    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_lw_in_diff   !< diffusion lw radiation
1043
1044!-- parameters required for RRTMG lower boundary condition
1045    REAL(wp)                   :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
1046    REAL(wp)                   :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
1047    REAL(wp)                   :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.
1048
1049!-- type for calculation of svf
1050    TYPE t_svf
1051        INTEGER(iwp)                               :: isurflt           !<
1052        INTEGER(iwp)                               :: isurfs            !<
1053        REAL(wp)                                   :: rsvf              !<
1054        REAL(wp)                                   :: rtransp           !<
1055    END TYPE
1056
1057!-- type for calculation of csf
1058    TYPE t_csf
1059        INTEGER(iwp)                               :: ip                !<
1060        INTEGER(iwp)                               :: itx               !<
1061        INTEGER(iwp)                               :: ity               !<
1062        INTEGER(iwp)                               :: itz               !<
1063        INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
1064        REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
1065                                                                        !< canopy sink factor for sky (-1)
1066    END TYPE
1067
1068!-- arrays storing the values of USM
1069    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
1070    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
1071    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
1072    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
1073
1074    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvf            !< array of sky view factor for each local surface
1075    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvft           !< array of sky view factor including transparency for each local surface
1076    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitrans         !< dsidir[isvfl,i] = path transmittance of i-th
1077                                                                        !< direction of direct solar irradiance per target surface
1078    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitransc        !< dtto per plant canopy box
1079    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsidir           !< dsidir[:,i] = unit vector of i-th
1080                                                                        !< direction of direct solar irradiance
1081    INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
1082    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
1083
1084    INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
1085    INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
1086    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
1087    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
1088    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
1089    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
1090    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
1091    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
1092    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
1093    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
1094    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
1095    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
1096    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
1097    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
1098    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box
1099
1100    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
1101    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw         !< array of lw radiation falling to local surface including radiation from reflections
1102    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir      !< array of direct sw radiation falling to local surface
1103    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif      !< array of diffuse sw radiation from sky and model boundary falling to local surface
1104    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif      !< array of diffuse lw radiation from sky and model boundary falling to local surface
1105   
1106                                                                        !< Outward radiation is only valid for nonvirtual surfaces
1107    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsl        !< array of reflected sw radiation for local surface in i-th reflection
1108    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutll        !< array of reflected + emitted lw radiation for local surface in i-th reflection
1109    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfouts         !< array of reflected sw radiation for all surfaces in i-th reflection
1110    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutl         !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
1111    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlg         !< global array of incoming lw radiation from plant canopy
1112    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw        !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1113    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1114    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
1115
1116!-- block variables needed for calculation of the plant canopy model inside the urban surface model
1117    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
1118    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
1119    REAL(wp), DIMENSION(:,:,:), POINTER            ::  sub_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
1120    REAL(wp), DIMENSION(:), POINTER                ::  sub_lad_g        !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
1121    REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
1122    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
1123    INTEGER(iwp)                                   ::  plantt_max
1124
1125!-- arrays and variables for calculation of svf and csf
1126    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
1127    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
1128    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
1129    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
1130    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
1131    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
1132    INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
1133    INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
1134    INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
1135    INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
1136    INTEGER(iwp), PARAMETER                        ::  gasize = 100000_iwp  !< initial size of growing arrays
1137    REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
1138    INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
1139    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
1140                                                                        !< needed only during calc_svf but must be here because it is
1141                                                                        !< shared between subroutines calc_svf and raytrace
1142    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
1143    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
1144    INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
1145
1146!-- temporary arrays for calculation of csf in raytracing
1147    INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
1148    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
1149    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
1150    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored
1151#if defined( __parallel )
1152    INTEGER(kind=MPI_ADDRESS_KIND), &
1153                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
1154    INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
1155    INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
1156#endif
1157    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
1158    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
1159    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
1160    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
1161    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_track_dist
1162    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
1163
1164!-- arrays for time averages
1165    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
1166    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
1167    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
1168    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
1169    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
1170    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
1171    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
1172    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
1173    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1174    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1175    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
1176    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
1177    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
1178    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
1179    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
1180    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
1181    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
1182
1183
1184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1185!-- Energy balance variables
1186!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1187!-- parameters of the land, roof and wall surfaces
1188    REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
1189    REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
1190
1191
1192    INTERFACE radiation_check_data_output
1193       MODULE PROCEDURE radiation_check_data_output
1194    END INTERFACE radiation_check_data_output
1195
1196    INTERFACE radiation_check_data_output_ts
1197       MODULE PROCEDURE radiation_check_data_output_ts
1198    END INTERFACE radiation_check_data_output_ts
1199
1200    INTERFACE radiation_check_data_output_pr
1201       MODULE PROCEDURE radiation_check_data_output_pr
1202    END INTERFACE radiation_check_data_output_pr
1203 
1204    INTERFACE radiation_check_parameters
1205       MODULE PROCEDURE radiation_check_parameters
1206    END INTERFACE radiation_check_parameters
1207 
1208    INTERFACE radiation_clearsky
1209       MODULE PROCEDURE radiation_clearsky
1210    END INTERFACE radiation_clearsky
1211 
1212    INTERFACE radiation_constant
1213       MODULE PROCEDURE radiation_constant
1214    END INTERFACE radiation_constant
1215 
1216    INTERFACE radiation_control
1217       MODULE PROCEDURE radiation_control
1218    END INTERFACE radiation_control
1219
1220    INTERFACE radiation_3d_data_averaging
1221       MODULE PROCEDURE radiation_3d_data_averaging
1222    END INTERFACE radiation_3d_data_averaging
1223
1224    INTERFACE radiation_data_output_2d
1225       MODULE PROCEDURE radiation_data_output_2d
1226    END INTERFACE radiation_data_output_2d
1227
1228    INTERFACE radiation_data_output_3d
1229       MODULE PROCEDURE radiation_data_output_3d
1230    END INTERFACE radiation_data_output_3d
1231
1232    INTERFACE radiation_data_output_mask
1233       MODULE PROCEDURE radiation_data_output_mask
1234    END INTERFACE radiation_data_output_mask
1235
1236    INTERFACE radiation_define_netcdf_grid
1237       MODULE PROCEDURE radiation_define_netcdf_grid
1238    END INTERFACE radiation_define_netcdf_grid
1239
1240    INTERFACE radiation_header
1241       MODULE PROCEDURE radiation_header
1242    END INTERFACE radiation_header
1243 
1244    INTERFACE radiation_init
1245       MODULE PROCEDURE radiation_init
1246    END INTERFACE radiation_init
1247
1248    INTERFACE radiation_parin
1249       MODULE PROCEDURE radiation_parin
1250    END INTERFACE radiation_parin
1251   
1252    INTERFACE radiation_rrtmg
1253       MODULE PROCEDURE radiation_rrtmg
1254    END INTERFACE radiation_rrtmg
1255
1256#if defined( __rrtmg )
1257    INTERFACE radiation_tendency
1258       MODULE PROCEDURE radiation_tendency
1259       MODULE PROCEDURE radiation_tendency_ij
1260    END INTERFACE radiation_tendency
1261#endif
1262
1263    INTERFACE radiation_rrd_local
1264       MODULE PROCEDURE radiation_rrd_local
1265    END INTERFACE radiation_rrd_local
1266
1267    INTERFACE radiation_wrd_local
1268       MODULE PROCEDURE radiation_wrd_local
1269    END INTERFACE radiation_wrd_local
1270
1271    INTERFACE radiation_interaction
1272       MODULE PROCEDURE radiation_interaction
1273    END INTERFACE radiation_interaction
1274
1275    INTERFACE radiation_interaction_init
1276       MODULE PROCEDURE radiation_interaction_init
1277    END INTERFACE radiation_interaction_init
1278 
1279    INTERFACE radiation_presimulate_solar_pos
1280       MODULE PROCEDURE radiation_presimulate_solar_pos
1281    END INTERFACE radiation_presimulate_solar_pos
1282
1283    INTERFACE radiation_calc_svf
1284       MODULE PROCEDURE radiation_calc_svf
1285    END INTERFACE radiation_calc_svf
1286
1287    INTERFACE radiation_write_svf
1288       MODULE PROCEDURE radiation_write_svf
1289    END INTERFACE radiation_write_svf
1290
1291    INTERFACE radiation_read_svf
1292       MODULE PROCEDURE radiation_read_svf
1293    END INTERFACE radiation_read_svf
1294
1295
1296    SAVE
1297
1298    PRIVATE
1299
1300!
1301!-- Public functions / NEEDS SORTING
1302    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
1303           radiation_check_data_output_ts,                                     &
1304           radiation_check_parameters, radiation_control,                      &
1305           radiation_header, radiation_init, radiation_parin,                  &
1306           radiation_3d_data_averaging,                                        &
1307           radiation_data_output_2d, radiation_data_output_3d,                 &
1308           radiation_define_netcdf_grid, radiation_wrd_local,                  &
1309           radiation_rrd_local, radiation_data_output_mask,                    &
1310           radiation_calc_svf, radiation_write_svf,                            &
1311           radiation_interaction, radiation_interaction_init,                  &
1312           radiation_read_svf, radiation_presimulate_solar_pos
1313
1314   
1315!
1316!-- Public variables and constants / NEEDS SORTING
1317    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
1318           emissivity, force_radiation_call, lat, lon, mrt_geom_human,         &
1319           mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
1320           rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
1321           rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
1322           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
1323           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
1324           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,           &
1325           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
1326           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
1327           idir, jdir, kdir, id, iz, iy, ix,                                   &
1328           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
1329           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
1330           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
1331           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
1332           radiation_interactions, startwall, startland, endland, endwall,     &
1333           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
1334           rad_sw_in_diff, rad_sw_in_dir
1335
1336
1337#if defined ( __rrtmg )
1338    PUBLIC radiation_tendency, rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
1339#endif
1340
1341 CONTAINS
1342
1343
1344!------------------------------------------------------------------------------!
1345! Description:
1346! ------------
1347!> This subroutine controls the calls of the radiation schemes
1348!------------------------------------------------------------------------------!
1349    SUBROUTINE radiation_control
1350 
1351 
1352       IMPLICIT NONE
1353
1354
1355       IF ( debug_output )  CALL debug_message( 'radiation_control', 'start' )
1356
1357
1358       SELECT CASE ( TRIM( radiation_scheme ) )
1359
1360          CASE ( 'constant' )
1361             CALL radiation_constant
1362         
1363          CASE ( 'clear-sky' ) 
1364             CALL radiation_clearsky
1365       
1366          CASE ( 'rrtmg' )
1367             CALL radiation_rrtmg
1368
1369          CASE DEFAULT
1370
1371       END SELECT
1372
1373       IF ( debug_output )  CALL debug_message( 'radiation_control', 'end' )
1374
1375    END SUBROUTINE radiation_control
1376
1377!------------------------------------------------------------------------------!
1378! Description:
1379! ------------
1380!> Check data output for radiation model
1381!------------------------------------------------------------------------------!
1382    SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
1383 
1384 
1385       USE control_parameters,                                                 &
1386           ONLY: data_output, message_string
1387
1388       IMPLICIT NONE
1389
1390       CHARACTER (LEN=*) ::  unit          !<
1391       CHARACTER (LEN=*) ::  variable      !<
1392
1393       INTEGER(iwp) :: i, k
1394       INTEGER(iwp) :: ilen
1395       CHARACTER(LEN=varnamelength) :: var          !< TRIM(variable)
1396
1397       var = TRIM(variable)
1398
1399!--    first process diractional variables
1400       IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
1401            var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
1402            var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
1403            var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
1404            var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
1405            var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
1406          IF ( .NOT.  radiation ) THEN
1407                message_string = 'output of "' // TRIM( var ) // '" require'&
1408                                 // 's radiation = .TRUE.'
1409                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1410          ENDIF
1411          unit = 'W/m2'
1412       ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
1413                 var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft'  .OR.             &
1414                 var(1:12) == 'rtm_surfalb_'  .OR.  var(1:13) == 'rtm_surfemis_'  ) THEN
1415          IF ( .NOT.  radiation ) THEN
1416                message_string = 'output of "' // TRIM( var ) // '" require'&
1417                                 // 's radiation = .TRUE.'
1418                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1419          ENDIF
1420          unit = '1'
1421       ELSE
1422!--       non-directional variables
1423          SELECT CASE ( TRIM( var ) )
1424             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
1425                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
1426                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
1427                   message_string = '"output of "' // TRIM( var ) // '" requi' // &
1428                                    'res radiation = .TRUE. and ' //              &
1429                                    'radiation_scheme = "rrtmg"'
1430                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
1431                ENDIF
1432                unit = 'K/h'
1433
1434             CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
1435                    'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
1436                    'rad_sw_out*')
1437                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
1438                   ! Workaround for masked output (calls with i=ilen=k=0)
1439                   unit = 'illegal'
1440                   RETURN
1441                ENDIF
1442                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
1443                   message_string = 'illegal value for data_output: "' //         &
1444                                    TRIM( var ) // '" & only 2d-horizontal ' //   &
1445                                    'cross sections are allowed for this value'
1446                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
1447                ENDIF
1448                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
1449                   IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
1450                        TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
1451                        TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
1452                        TRIM( var ) == 'rrtm_asdir*'      )                       &
1453                   THEN
1454                      message_string = 'output of "' // TRIM( var ) // '" require'&
1455                                       // 's radiation = .TRUE. and radiation_sch'&
1456                                       // 'eme = "rrtmg"'
1457                      CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
1458                   ENDIF
1459                ENDIF
1460
1461                IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
1462                IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
1463                IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
1464                IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
1465                IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
1466                IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
1467                IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''
1468                IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
1469                IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
1470                IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
1471
1472             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
1473                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
1474                IF ( .NOT.  radiation ) THEN
1475                   message_string = 'output of "' // TRIM( var ) // '" require'&
1476                                    // 's radiation = .TRUE.'
1477                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1478                ENDIF
1479                unit = 'W'
1480
1481             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
1482                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
1483                   ! Workaround for masked output (calls with i=ilen=k=0)
1484                   unit = 'illegal'
1485                   RETURN
1486                ENDIF
1487
1488                IF ( .NOT.  radiation ) THEN
1489                   message_string = 'output of "' // TRIM( var ) // '" require'&
1490                                    // 's radiation = .TRUE.'
1491                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1492                ENDIF
1493                IF ( mrt_nlevels == 0 ) THEN
1494                   message_string = 'output of "' // TRIM( var ) // '" require'&
1495                                    // 's mrt_nlevels > 0'
1496                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
1497                ENDIF
1498                IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
1499                   message_string = 'output of "' // TRIM( var ) // '" require'&
1500                                    // 's rtm_mrt_sw = .TRUE.'
1501                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
1502                ENDIF
1503                IF ( TRIM( var ) == 'rtm_mrt' ) THEN
1504                   unit = 'K'
1505                ELSE
1506                   unit = 'W m-2'
1507                ENDIF
1508
1509             CASE DEFAULT
1510                unit = 'illegal'
1511
1512          END SELECT
1513       ENDIF
1514
1515    END SUBROUTINE radiation_check_data_output
1516
1517
1518!------------------------------------------------------------------------------!
1519! Description:
1520! ------------
1521!> Set module-specific timeseries units and labels
1522!------------------------------------------------------------------------------!
1523 SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )
1524
1525
1526    INTEGER(iwp),      INTENT(IN)     ::  dots_max
1527    INTEGER(iwp),      INTENT(INOUT)  ::  dots_num
1528
1529!
1530!-- Next line is just to avoid compiler warning about unused variable.
1531    IF ( dots_max == 0 )  CONTINUE
1532
1533!
1534!-- Temporary solution to add LSM and radiation time series to the default
1535!-- output
1536    IF ( land_surface  .OR.  radiation )  THEN
1537       IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
1538          dots_num = dots_num + 15
1539       ELSE
1540          dots_num = dots_num + 11
1541       ENDIF
1542    ENDIF
1543
1544
1545 END SUBROUTINE radiation_check_data_output_ts
1546
1547!------------------------------------------------------------------------------!
1548! Description:
1549! ------------
1550!> Check data output of profiles for radiation model
1551!------------------------------------------------------------------------------! 
1552    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
1553               dopr_unit )
1554 
1555       USE arrays_3d,                                                          &
1556           ONLY: zu
1557
1558       USE control_parameters,                                                 &
1559           ONLY: data_output_pr, message_string
1560
1561       USE indices
1562
1563       USE profil_parameter
1564
1565       USE statistics
1566
1567       IMPLICIT NONE
1568   
1569       CHARACTER (LEN=*) ::  unit      !<
1570       CHARACTER (LEN=*) ::  variable  !<
1571       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
1572 
1573       INTEGER(iwp) ::  var_count     !<
1574
1575       SELECT CASE ( TRIM( variable ) )
1576       
1577         CASE ( 'rad_net' )
1578             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
1579             THEN
1580                message_string = 'data_output_pr = ' //                        &
1581                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1582                                 'not available for radiation = .FALSE. or ' //&
1583                                 'radiation_scheme = "constant"'
1584                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1585             ELSE
1586                dopr_index(var_count) = 99
1587                dopr_unit  = 'W/m2'
1588                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
1589                unit = dopr_unit
1590             ENDIF
1591
1592          CASE ( 'rad_lw_in' )
1593             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
1594             THEN
1595                message_string = 'data_output_pr = ' //                        &
1596                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1597                                 'not available for radiation = .FALSE. or ' //&
1598                                 'radiation_scheme = "constant"'
1599                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1600             ELSE
1601                dopr_index(var_count) = 100
1602                dopr_unit  = 'W/m2'
1603                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
1604                unit = dopr_unit
1605             ENDIF
1606
1607          CASE ( 'rad_lw_out' )
1608             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
1609             THEN
1610                message_string = 'data_output_pr = ' //                        &
1611                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1612                                 'not available for radiation = .FALSE. or ' //&
1613                                 'radiation_scheme = "constant"'
1614                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1615             ELSE
1616                dopr_index(var_count) = 101
1617                dopr_unit  = 'W/m2'
1618                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
1619                unit = dopr_unit   
1620             ENDIF
1621
1622          CASE ( 'rad_sw_in' )
1623             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
1624             THEN
1625                message_string = 'data_output_pr = ' //                        &
1626                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1627                                 'not available for radiation = .FALSE. or ' //&
1628                                 'radiation_scheme = "constant"'
1629                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1630             ELSE
1631                dopr_index(var_count) = 102
1632                dopr_unit  = 'W/m2'
1633                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
1634                unit = dopr_unit
1635             ENDIF
1636
1637          CASE ( 'rad_sw_out')
1638             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
1639             THEN
1640                message_string = 'data_output_pr = ' //                        &
1641                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1642                                 'not available for radiation = .FALSE. or ' //&
1643                                 'radiation_scheme = "constant"'
1644                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1645             ELSE
1646                dopr_index(var_count) = 103
1647                dopr_unit  = 'W/m2'
1648                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
1649                unit = dopr_unit
1650             ENDIF
1651
1652          CASE ( 'rad_lw_cs_hr' )
1653             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
1654             THEN
1655                message_string = 'data_output_pr = ' //                        &
1656                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1657                                 'not available for radiation = .FALSE. or ' //&
1658                                 'radiation_scheme /= "rrtmg"'
1659                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1660             ELSE
1661                dopr_index(var_count) = 104
1662                dopr_unit  = 'K/h'
1663                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
1664                unit = dopr_unit
1665             ENDIF
1666
1667          CASE ( 'rad_lw_hr' )
1668             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
1669             THEN
1670                message_string = 'data_output_pr = ' //                        &
1671                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1672                                 'not available for radiation = .FALSE. or ' //&
1673                                 'radiation_scheme /= "rrtmg"'
1674                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1675             ELSE
1676                dopr_index(var_count) = 105
1677                dopr_unit  = 'K/h'
1678                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
1679                unit = dopr_unit
1680             ENDIF
1681
1682          CASE ( 'rad_sw_cs_hr' )
1683             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
1684             THEN
1685                message_string = 'data_output_pr = ' //                        &
1686                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1687                                 'not available for radiation = .FALSE. or ' //&
1688                                 'radiation_scheme /= "rrtmg"'
1689                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1690             ELSE
1691                dopr_index(var_count) = 106
1692                dopr_unit  = 'K/h'
1693                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
1694                unit = dopr_unit
1695             ENDIF
1696
1697          CASE ( 'rad_sw_hr' )
1698             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
1699             THEN
1700                message_string = 'data_output_pr = ' //                        &
1701                                 TRIM( data_output_pr(var_count) ) // ' is' // &
1702                                 'not available for radiation = .FALSE. or ' //&
1703                                 'radiation_scheme /= "rrtmg"'
1704                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1705             ELSE
1706                dopr_index(var_count) = 107
1707                dopr_unit  = 'K/h'
1708                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
1709                unit = dopr_unit
1710             ENDIF
1711
1712
1713          CASE DEFAULT
1714             unit = 'illegal'
1715
1716       END SELECT
1717
1718
1719    END SUBROUTINE radiation_check_data_output_pr
1720 
1721 
1722!------------------------------------------------------------------------------!
1723! Description:
1724! ------------
1725!> Check parameters routine for radiation model
1726!------------------------------------------------------------------------------!
1727    SUBROUTINE radiation_check_parameters
1728
1729       USE control_parameters,                                                 &
1730           ONLY: land_surface, message_string, urban_surface
1731
1732       USE netcdf_data_input_mod,                                              &
1733           ONLY:  input_pids_static                 
1734   
1735       IMPLICIT NONE
1736       
1737!
1738!--    In case no urban-surface or land-surface model is applied, usage of
1739!--    a radiation model make no sense.         
1740       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
1741          message_string = 'Usage of radiation module is only allowed if ' //  &
1742                           'land-surface and/or urban-surface model is applied.'
1743          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1744       ENDIF
1745
1746       IF ( radiation_scheme /= 'constant'   .AND.                             &
1747            radiation_scheme /= 'clear-sky'  .AND.                             &
1748            radiation_scheme /= 'rrtmg' )  THEN
1749          message_string = 'unknown radiation_scheme = '//                     &
1750                           TRIM( radiation_scheme )
1751          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
1752       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
1753#if ! defined ( __rrtmg )
1754          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
1755                           'compilation of PALM with pre-processor ' //        &
1756                           'directive -D__rrtmg'
1757          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
1758#endif
1759#if defined ( __rrtmg ) && ! defined( __netcdf )
1760          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
1761                           'the use of NetCDF (preprocessor directive ' //     &
1762                           '-D__netcdf'
1763          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
1764#endif
1765
1766       ENDIF
1767!
1768!--    Checks performed only if data is given via namelist only.
1769       IF ( .NOT. input_pids_static )  THEN
1770          IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.          &
1771               radiation_scheme == 'clear-sky')  THEN
1772             message_string = 'radiation_scheme = "clear-sky" in combination'//& 
1773                              'with albedo_type = 0 requires setting of'//     &
1774                              'albedo /= 9999999.9'
1775             CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
1776          ENDIF
1777
1778          IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.     &
1779             ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
1780          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
1781             ) ) THEN
1782             message_string = 'radiation_scheme = "rrtmg" in combination' //   & 
1783                              'with albedo_type = 0 requires setting of ' //   &
1784                              'albedo_lw_dif /= 9999999.9' //                  &
1785                              'albedo_lw_dir /= 9999999.9' //                  &
1786                              'albedo_sw_dif /= 9999999.9 and' //              &
1787                              'albedo_sw_dir /= 9999999.9'
1788             CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
1789          ENDIF
1790       ENDIF
1791!
1792!--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
1793#if defined( __parallel )     
1794       IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
1795          message_string = 'rad_angular_discretization can only be used ' //  &
1796                           'together with raytrace_mpi_rma or when ' //  &
1797                           'no parallelization is applied.'
1798          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1799       ENDIF
1800#endif
1801
1802       IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
1803            average_radiation ) THEN
1804          message_string = 'average_radiation = .T. with radiation_scheme'//   & 
1805                           '= "rrtmg" in combination cloud_droplets = .T.'//   &
1806                           'is not implementd'
1807          CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
1808       ENDIF
1809
1810!
1811!--    Incialize svf normalization reporting histogram
1812       svfnorm_report_num = 1
1813       DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp          &
1814                   .AND. svfnorm_report_num <= 30 )
1815          svfnorm_report_num = svfnorm_report_num + 1
1816       ENDDO
1817       svfnorm_report_num = svfnorm_report_num - 1
1818!
1819!--    Check for dt_radiation
1820       IF ( dt_radiation <= 0.0 )  THEN
1821          message_string = 'dt_radiation must be > 0.0' 
1822          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 
1823       ENDIF
1824 
1825    END SUBROUTINE radiation_check_parameters
1826 
1827 
1828!------------------------------------------------------------------------------!
1829! Description:
1830! ------------
1831!> Initialization of the radiation model
1832!------------------------------------------------------------------------------!
1833    SUBROUTINE radiation_init
1834   
1835       IMPLICIT NONE
1836
1837       INTEGER(iwp) ::  i         !< running index x-direction
1838       INTEGER(iwp) ::  ioff      !< offset in x between surface element reference grid point in atmosphere and actual surface
1839       INTEGER(iwp) ::  j         !< running index y-direction
1840       INTEGER(iwp) ::  joff      !< offset in y between surface element reference grid point in atmosphere and actual surface
1841       INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
1842       INTEGER(iwp) ::  m         !< running index for surface elements
1843#if defined( __rrtmg )
1844       INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
1845#endif
1846
1847
1848       IF ( debug_output )  CALL debug_message( 'radiation_init', 'start' )
1849!
1850!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
1851!--    The namelist parameter radiation_interactions_on can override this behavior.
1852!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
1853!--    init_surface_arrays.)
1854       IF ( radiation_interactions_on )  THEN
1855          IF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
1856             radiation_interactions    = .TRUE.
1857             average_radiation         = .TRUE.
1858          ELSE
1859             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
1860                                                   !< calculations necessary in case of flat surface
1861          ENDIF
1862       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
1863          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
1864                           'vertical surfaces and/or trees exist. The model will run ' // &
1865                           'without RTM (no shadows, no radiation reflections)'
1866          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
1867       ENDIF
1868!
1869!--    If required, initialize radiation interactions between surfaces
1870!--    via sky-view factors. This must be done before radiation is initialized.
1871       IF ( radiation_interactions )  CALL radiation_interaction_init
1872!
1873!--    Allocate array for storing the surface net radiation
1874       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.                      &
1875                  surf_lsm_h%ns > 0  )   THEN
1876          ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
1877          surf_lsm_h%rad_net = 0.0_wp 
1878       ENDIF
1879       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
1880                  surf_usm_h%ns > 0  )  THEN
1881          ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
1882          surf_usm_h%rad_net = 0.0_wp 
1883       ENDIF
1884       DO  l = 0, 3
1885          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.                &
1886                     surf_lsm_v(l)%ns > 0  )  THEN
1887             ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
1888             surf_lsm_v(l)%rad_net = 0.0_wp 
1889          ENDIF
1890          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.                &
1891                     surf_usm_v(l)%ns > 0  )  THEN
1892             ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
1893             surf_usm_v(l)%rad_net = 0.0_wp 
1894          ENDIF
1895       ENDDO
1896
1897
1898!
1899!--    Allocate array for storing the surface longwave (out) radiation change
1900       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.          &
1901                  surf_lsm_h%ns > 0  )   THEN
1902          ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
1903          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp 
1904       ENDIF
1905       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
1906                  surf_usm_h%ns > 0  )  THEN
1907          ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
1908          surf_usm_h%rad_lw_out_change_0 = 0.0_wp 
1909       ENDIF
1910       DO  l = 0, 3
1911          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.    &
1912                     surf_lsm_v(l)%ns > 0  )  THEN
1913             ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
1914             surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp 
1915          ENDIF
1916          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.    &
1917                     surf_usm_v(l)%ns > 0  )  THEN
1918             ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
1919             surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp 
1920          ENDIF
1921       ENDDO
1922
1923!
1924!--    Allocate surface arrays for incoming/outgoing short/longwave radiation
1925       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.                    &
1926                  surf_lsm_h%ns > 0  )   THEN
1927          ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
1928          ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
1929          ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
1930          ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
1931          ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
1932          ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
1933          ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
1934          ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
1935          ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
1936          ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
1937          ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
1938          surf_lsm_h%rad_sw_in  = 0.0_wp 
1939          surf_lsm_h%rad_sw_out = 0.0_wp 
1940          surf_lsm_h%rad_sw_dir = 0.0_wp 
1941          surf_lsm_h%rad_sw_dif = 0.0_wp 
1942          surf_lsm_h%rad_sw_ref = 0.0_wp 
1943          surf_lsm_h%rad_sw_res = 0.0_wp 
1944          surf_lsm_h%rad_lw_in  = 0.0_wp 
1945          surf_lsm_h%rad_lw_out = 0.0_wp 
1946          surf_lsm_h%rad_lw_dif = 0.0_wp 
1947          surf_lsm_h%rad_lw_ref = 0.0_wp 
1948          surf_lsm_h%rad_lw_res = 0.0_wp 
1949       ENDIF
1950       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
1951                  surf_usm_h%ns > 0  )  THEN
1952          ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
1953          ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
1954          ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
1955          ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
1956          ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
1957          ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
1958          ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
1959          ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
1960          ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
1961          ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
1962          ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
1963          surf_usm_h%rad_sw_in  = 0.0_wp 
1964          surf_usm_h%rad_sw_out = 0.0_wp 
1965          surf_usm_h%rad_sw_dir = 0.0_wp 
1966          surf_usm_h%rad_sw_dif = 0.0_wp 
1967          surf_usm_h%rad_sw_ref = 0.0_wp 
1968          surf_usm_h%rad_sw_res = 0.0_wp 
1969          surf_usm_h%rad_lw_in  = 0.0_wp 
1970          surf_usm_h%rad_lw_out = 0.0_wp 
1971          surf_usm_h%rad_lw_dif = 0.0_wp 
1972          surf_usm_h%rad_lw_ref = 0.0_wp 
1973          surf_usm_h%rad_lw_res = 0.0_wp 
1974       ENDIF
1975       DO  l = 0, 3
1976          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.              &
1977                     surf_lsm_v(l)%ns > 0  )  THEN
1978             ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
1979             ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
1980             ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
1981             ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
1982             ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
1983             ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
1984
1985             ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
1986             ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
1987             ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
1988             ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
1989             ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
1990
1991             surf_lsm_v(l)%rad_sw_in  = 0.0_wp 
1992             surf_lsm_v(l)%rad_sw_out = 0.0_wp
1993             surf_lsm_v(l)%rad_sw_dir = 0.0_wp
1994             surf_lsm_v(l)%rad_sw_dif = 0.0_wp
1995             surf_lsm_v(l)%rad_sw_ref = 0.0_wp
1996             surf_lsm_v(l)%rad_sw_res = 0.0_wp
1997
1998             surf_lsm_v(l)%rad_lw_in  = 0.0_wp 
1999             surf_lsm_v(l)%rad_lw_out = 0.0_wp 
2000             surf_lsm_v(l)%rad_lw_dif = 0.0_wp 
2001             surf_lsm_v(l)%rad_lw_ref = 0.0_wp 
2002             surf_lsm_v(l)%rad_lw_res = 0.0_wp 
2003          ENDIF
2004          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.              &
2005                     surf_usm_v(l)%ns > 0  )  THEN
2006             ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
2007             ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
2008             ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
2009             ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
2010             ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
2011             ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
2012             ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
2013             ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
2014             ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
2015             ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
2016             ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
2017             surf_usm_v(l)%rad_sw_in  = 0.0_wp 
2018             surf_usm_v(l)%rad_sw_out = 0.0_wp
2019             surf_usm_v(l)%rad_sw_dir = 0.0_wp
2020             surf_usm_v(l)%rad_sw_dif = 0.0_wp
2021             surf_usm_v(l)%rad_sw_ref = 0.0_wp
2022             surf_usm_v(l)%rad_sw_res = 0.0_wp
2023             surf_usm_v(l)%rad_lw_in  = 0.0_wp 
2024             surf_usm_v(l)%rad_lw_out = 0.0_wp 
2025             surf_usm_v(l)%rad_lw_dif = 0.0_wp 
2026             surf_usm_v(l)%rad_lw_ref = 0.0_wp 
2027             surf_usm_v(l)%rad_lw_res = 0.0_wp 
2028          ENDIF
2029       ENDDO
2030!
2031!--    Fix net radiation in case of radiation_scheme = 'constant'
2032       IF ( radiation_scheme == 'constant' )  THEN
2033          IF ( ALLOCATED( surf_lsm_h%rad_net ) )                               &
2034             surf_lsm_h%rad_net    = net_radiation
2035          IF ( ALLOCATED( surf_usm_h%rad_net ) )                               &
2036             surf_usm_h%rad_net    = net_radiation
2037!
2038!--       Todo: weight with inclination angle
2039          DO  l = 0, 3
2040             IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )                         &
2041                surf_lsm_v(l)%rad_net = net_radiation
2042             IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )                         &
2043                surf_usm_v(l)%rad_net = net_radiation
2044          ENDDO
2045!          radiation = .FALSE.
2046!
2047!--    Calculate orbital constants
2048       ELSE
2049          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
2050          decl_2 = 2.0_wp * pi / 365.0_wp
2051          decl_3 = decl_2 * 81.0_wp
2052          lat    = latitude * pi / 180.0_wp
2053          lon    = longitude * pi / 180.0_wp
2054       ENDIF
2055
2056       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
2057            radiation_scheme == 'constant')  THEN
2058
2059
2060!
2061!--       Allocate arrays for incoming/outgoing short/longwave radiation
2062          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
2063             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
2064          ENDIF
2065          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
2066             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
2067          ENDIF
2068
2069          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
2070             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
2071          ENDIF
2072          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
2073             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
2074          ENDIF
2075
2076!
2077!--       Allocate average arrays for incoming/outgoing short/longwave radiation
2078          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
2079             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
2080          ENDIF
2081          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
2082             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
2083          ENDIF
2084
2085          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
2086             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
2087          ENDIF
2088          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
2089             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
2090          ENDIF
2091!
2092!--       Allocate arrays for broadband albedo, and level 1 initialization
2093!--       via namelist paramter, unless not already allocated.
2094          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
2095             ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns)     )
2096             surf_lsm_h%albedo    = albedo
2097          ENDIF
2098          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
2099             ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)     )
2100             surf_usm_h%albedo    = albedo
2101          ENDIF
2102
2103          DO  l = 0, 3
2104             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
2105                ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
2106                surf_lsm_v(l)%albedo = albedo
2107             ENDIF
2108             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
2109                ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
2110                surf_usm_v(l)%albedo = albedo
2111             ENDIF
2112          ENDDO
2113!
2114!--       Level 2 initialization of broadband albedo via given albedo_type.
2115!--       Only if albedo_type is non-zero. In case of urban surface and
2116!--       input data is read from ASCII file, albedo_type will be zero, so that
2117!--       albedo won't be overwritten.
2118          DO  m = 1, surf_lsm_h%ns
2119             IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) /= 0 )                &
2120                surf_lsm_h%albedo(ind_veg_wall,m) =                            &
2121                           albedo_pars(2,surf_lsm_h%albedo_type(ind_veg_wall,m))
2122             IF ( surf_lsm_h%albedo_type(ind_pav_green,m) /= 0 )               &
2123                surf_lsm_h%albedo(ind_pav_green,m) =                           &
2124                           albedo_pars(2,surf_lsm_h%albedo_type(ind_pav_green,m))
2125             IF ( surf_lsm_h%albedo_type(ind_wat_win,m) /= 0 )                 &
2126                surf_lsm_h%albedo(ind_wat_win,m) =                             &
2127                           albedo_pars(2,surf_lsm_h%albedo_type(ind_wat_win,m))
2128          ENDDO
2129          DO  m = 1, surf_usm_h%ns
2130             IF ( surf_usm_h%albedo_type(ind_veg_wall,m) /= 0 )                &
2131                surf_usm_h%albedo(ind_veg_wall,m) =                            &
2132                           albedo_pars(2,surf_usm_h%albedo_type(ind_veg_wall,m))
2133             IF ( surf_usm_h%albedo_type(ind_pav_green,m) /= 0 )               &
2134                surf_usm_h%albedo(ind_pav_green,m) =                           &
2135                           albedo_pars(2,surf_usm_h%albedo_type(ind_pav_green,m))
2136             IF ( surf_usm_h%albedo_type(ind_wat_win,m) /= 0 )                 &
2137                surf_usm_h%albedo(ind_wat_win,m) =                             &
2138                           albedo_pars(2,surf_usm_h%albedo_type(ind_wat_win,m))
2139          ENDDO
2140
2141          DO  l = 0, 3
2142             DO  m = 1, surf_lsm_v(l)%ns
2143                IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) /= 0 )          &
2144                   surf_lsm_v(l)%albedo(ind_veg_wall,m) =                      &
2145                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_veg_wall,m))
2146                IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) /= 0 )         &
2147                   surf_lsm_v(l)%albedo(ind_pav_green,m) =                     &
2148                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_pav_green,m))
2149                IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) /= 0 )           &
2150                   surf_lsm_v(l)%albedo(ind_wat_win,m) =                       &
2151                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_wat_win,m))
2152             ENDDO
2153             DO  m = 1, surf_usm_v(l)%ns
2154                IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) /= 0 )          &
2155                   surf_usm_v(l)%albedo(ind_veg_wall,m) =                      &
2156                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_veg_wall,m))
2157                IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) /= 0 )         &
2158                   surf_usm_v(l)%albedo(ind_pav_green,m) =                     &
2159                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_pav_green,m))
2160                IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) /= 0 )           &
2161                   surf_usm_v(l)%albedo(ind_wat_win,m) =                       &
2162                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_wat_win,m))
2163             ENDDO
2164          ENDDO
2165
2166!
2167!--       Level 3 initialization at grid points where albedo type is zero.
2168!--       This case, albedo is taken from file. In case of constant radiation
2169!--       or clear sky, only broadband albedo is given.
2170          IF ( albedo_pars_f%from_file )  THEN
2171!
2172!--          Horizontal surfaces
2173             DO  m = 1, surf_lsm_h%ns
2174                i = surf_lsm_h%i(m)
2175                j = surf_lsm_h%j(m)
2176                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
2177                   IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) == 0 )          &
2178                      surf_lsm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2179                   IF ( surf_lsm_h%albedo_type(ind_pav_green,m) == 0 )         &
2180                      surf_lsm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2181                   IF ( surf_lsm_h%albedo_type(ind_wat_win,m) == 0 )           &
2182                      surf_lsm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2183                ENDIF
2184             ENDDO
2185             DO  m = 1, surf_usm_h%ns
2186                i = surf_usm_h%i(m)
2187                j = surf_usm_h%j(m)
2188                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
2189                   IF ( surf_usm_h%albedo_type(ind_veg_wall,m) == 0 )          &
2190                      surf_usm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2191                   IF ( surf_usm_h%albedo_type(ind_pav_green,m) == 0 )         &
2192                      surf_usm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2193                   IF ( surf_usm_h%albedo_type(ind_wat_win,m) == 0 )           &
2194                      surf_usm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2195                ENDIF
2196             ENDDO
2197!
2198!--          Vertical surfaces           
2199             DO  l = 0, 3
2200
2201                ioff = surf_lsm_v(l)%ioff
2202                joff = surf_lsm_v(l)%joff
2203                DO  m = 1, surf_lsm_v(l)%ns
2204                   i = surf_lsm_v(l)%i(m) + ioff
2205                   j = surf_lsm_v(l)%j(m) + joff
2206                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
2207                      IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) == 0 )    &
2208                         surf_lsm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2209                      IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) == 0 )   &
2210                         surf_lsm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2211                      IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) == 0 )     &
2212                         surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2213                   ENDIF
2214                ENDDO
2215
2216                ioff = surf_usm_v(l)%ioff
2217                joff = surf_usm_v(l)%joff
2218                DO  m = 1, surf_usm_h%ns
2219                   i = surf_usm_h%i(m) + joff
2220                   j = surf_usm_h%j(m) + joff
2221                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
2222                      IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) == 0 )    &
2223                         surf_usm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2224                      IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) == 0 )   &
2225                         surf_usm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2226                      IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) == 0 )     &
2227                         surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2228                   ENDIF
2229                ENDDO
2230             ENDDO
2231
2232          ENDIF 
2233!
2234!--    Initialization actions for RRTMG
2235       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
2236#if defined ( __rrtmg )
2237!
2238!--       Allocate albedos for short/longwave radiation, horizontal surfaces
2239!--       for wall/green/window (USM) or vegetation/pavement/water surfaces
2240!--       (LSM).
2241          ALLOCATE ( surf_lsm_h%aldif(0:2,1:surf_lsm_h%ns)       )
2242          ALLOCATE ( surf_lsm_h%aldir(0:2,1:surf_lsm_h%ns)       )
2243          ALLOCATE ( surf_lsm_h%asdif(0:2,1:surf_lsm_h%ns)       )
2244          ALLOCATE ( surf_lsm_h%asdir(0:2,1:surf_lsm_h%ns)       )
2245          ALLOCATE ( surf_lsm_h%rrtm_aldif(0:2,1:surf_lsm_h%ns)  )
2246          ALLOCATE ( surf_lsm_h%rrtm_aldir(0:2,1:surf_lsm_h%ns)  )
2247          ALLOCATE ( surf_lsm_h%rrtm_asdif(0:2,1:surf_lsm_h%ns)  )
2248          ALLOCATE ( surf_lsm_h%rrtm_asdir(0:2,1:surf_lsm_h%ns)  )
2249
2250          ALLOCATE ( surf_usm_h%aldif(0:2,1:surf_usm_h%ns)       )
2251          ALLOCATE ( surf_usm_h%aldir(0:2,1:surf_usm_h%ns)       )
2252          ALLOCATE ( surf_usm_h%asdif(0:2,1:surf_usm_h%ns)       )
2253          ALLOCATE ( surf_usm_h%asdir(0:2,1:surf_usm_h%ns)       )
2254          ALLOCATE ( surf_usm_h%rrtm_aldif(0:2,1:surf_usm_h%ns)  )
2255          ALLOCATE ( surf_usm_h%rrtm_aldir(0:2,1:surf_usm_h%ns)  )
2256          ALLOCATE ( surf_usm_h%rrtm_asdif(0:2,1:surf_usm_h%ns)  )
2257          ALLOCATE ( surf_usm_h%rrtm_asdir(0:2,1:surf_usm_h%ns)  )
2258
2259!
2260!--       Allocate broadband albedo (temporary for the current radiation
2261!--       implementations)
2262          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
2263             ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns)     )
2264          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
2265             ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)     )
2266
2267!
2268!--       Allocate albedos for short/longwave radiation, vertical surfaces
2269          DO  l = 0, 3
2270
2271             ALLOCATE ( surf_lsm_v(l)%aldif(0:2,1:surf_lsm_v(l)%ns)      )
2272             ALLOCATE ( surf_lsm_v(l)%aldir(0:2,1:surf_lsm_v(l)%ns)      )
2273             ALLOCATE ( surf_lsm_v(l)%asdif(0:2,1:surf_lsm_v(l)%ns)      )
2274             ALLOCATE ( surf_lsm_v(l)%asdir(0:2,1:surf_lsm_v(l)%ns)      )
2275
2276             ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(0:2,1:surf_lsm_v(l)%ns) )
2277             ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(0:2,1:surf_lsm_v(l)%ns) )
2278             ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(0:2,1:surf_lsm_v(l)%ns) )
2279             ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(0:2,1:surf_lsm_v(l)%ns) )
2280
2281             ALLOCATE ( surf_usm_v(l)%aldif(0:2,1:surf_usm_v(l)%ns)      )
2282             ALLOCATE ( surf_usm_v(l)%aldir(0:2,1:surf_usm_v(l)%ns)      )
2283             ALLOCATE ( surf_usm_v(l)%asdif(0:2,1:surf_usm_v(l)%ns)      )
2284             ALLOCATE ( surf_usm_v(l)%asdir(0:2,1:surf_usm_v(l)%ns)      )
2285
2286             ALLOCATE ( surf_usm_v(l)%rrtm_aldif(0:2,1:surf_usm_v(l)%ns) )
2287             ALLOCATE ( surf_usm_v(l)%rrtm_aldir(0:2,1:surf_usm_v(l)%ns) )
2288             ALLOCATE ( surf_usm_v(l)%rrtm_asdif(0:2,1:surf_usm_v(l)%ns) )
2289             ALLOCATE ( surf_usm_v(l)%rrtm_asdir(0:2,1:surf_usm_v(l)%ns) )
2290!
2291!--          Allocate broadband albedo (temporary for the current radiation
2292!--          implementations)
2293             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                    &
2294                ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
2295             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                    &
2296                ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
2297
2298          ENDDO
2299!
2300!--       Level 1 initialization of spectral albedos via namelist
2301!--       paramters. Please note, this case all surface tiles are initialized
2302!--       the same.
2303          IF ( surf_lsm_h%ns > 0 )  THEN
2304             surf_lsm_h%aldif  = albedo_lw_dif
2305             surf_lsm_h%aldir  = albedo_lw_dir
2306             surf_lsm_h%asdif  = albedo_sw_dif
2307             surf_lsm_h%asdir  = albedo_sw_dir
2308             surf_lsm_h%albedo = albedo_sw_dif
2309          ENDIF
2310          IF ( surf_usm_h%ns > 0 )  THEN
2311             IF ( surf_usm_h%albedo_from_ascii )  THEN
2312                surf_usm_h%aldif  = surf_usm_h%albedo
2313                surf_usm_h%aldir  = surf_usm_h%albedo
2314                surf_usm_h%asdif  = surf_usm_h%albedo
2315                surf_usm_h%asdir  = surf_usm_h%albedo
2316             ELSE
2317                surf_usm_h%aldif  = albedo_lw_dif
2318                surf_usm_h%aldir  = albedo_lw_dir
2319                surf_usm_h%asdif  = albedo_sw_dif
2320                surf_usm_h%asdir  = albedo_sw_dir
2321                surf_usm_h%albedo = albedo_sw_dif
2322             ENDIF
2323          ENDIF
2324
2325          DO  l = 0, 3
2326
2327             IF ( surf_lsm_v(l)%ns > 0 )  THEN
2328                surf_lsm_v(l)%aldif  = albedo_lw_dif
2329                surf_lsm_v(l)%aldir  = albedo_lw_dir
2330                surf_lsm_v(l)%asdif  = albedo_sw_dif
2331                surf_lsm_v(l)%asdir  = albedo_sw_dir
2332                surf_lsm_v(l)%albedo = albedo_sw_dif
2333             ENDIF
2334
2335             IF ( surf_usm_v(l)%ns > 0 )  THEN
2336                IF ( surf_usm_v(l)%albedo_from_ascii )  THEN
2337                   surf_usm_v(l)%aldif  = surf_usm_v(l)%albedo
2338                   surf_usm_v(l)%aldir  = surf_usm_v(l)%albedo
2339                   surf_usm_v(l)%asdif  = surf_usm_v(l)%albedo
2340                   surf_usm_v(l)%asdir  = surf_usm_v(l)%albedo
2341                ELSE
2342                   surf_usm_v(l)%aldif  = albedo_lw_dif
2343                   surf_usm_v(l)%aldir  = albedo_lw_dir
2344                   surf_usm_v(l)%asdif  = albedo_sw_dif
2345                   surf_usm_v(l)%asdir  = albedo_sw_dir
2346                ENDIF
2347             ENDIF
2348          ENDDO
2349
2350!
2351!--       Level 2 initialization of spectral albedos via albedo_type.
2352!--       Please note, for natural- and urban-type surfaces, a tile approach
2353!--       is applied so that the resulting albedo is calculated via the weighted
2354!--       average of respective surface fractions.
2355          DO  m = 1, surf_lsm_h%ns
2356!
2357!--          Spectral albedos for vegetation/pavement/water surfaces
2358             DO  ind_type = 0, 2
2359                IF ( surf_lsm_h%albedo_type(ind_type,m) /= 0 )  THEN
2360                   surf_lsm_h%aldif(ind_type,m) =                              &
2361                               albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
2362                   surf_lsm_h%asdif(ind_type,m) =                              &
2363                               albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2364                   surf_lsm_h%aldir(ind_type,m) =                              &
2365                               albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
2366                   surf_lsm_h%asdir(ind_type,m) =                              &
2367                               albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2368                   surf_lsm_h%albedo(ind_type,m) =                             &
2369                               albedo_pars(2,surf_lsm_h%albedo_type(ind_type,m))
2370                ENDIF
2371             ENDDO
2372
2373          ENDDO
2374!
2375!--       For urban surface only if albedo has not been already initialized
2376!--       in the urban-surface model via the ASCII file.
2377          IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
2378             DO  m = 1, surf_usm_h%ns
2379!
2380!--             Spectral albedos for wall/green/window surfaces
2381                DO  ind_type = 0, 2
2382                   IF ( surf_usm_h%albedo_type(ind_type,m) /= 0 )  THEN
2383                      surf_usm_h%aldif(ind_type,m) =                           &
2384                               albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
2385                      surf_usm_h%asdif(ind_type,m) =                           &
2386                               albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2387                      surf_usm_h%aldir(ind_type,m) =                           &
2388                               albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
2389                      surf_usm_h%asdir(ind_type,m) =                           &
2390                               albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2391                      surf_usm_h%albedo(ind_type,m) =                          &
2392                               albedo_pars(2,surf_usm_h%albedo_type(ind_type,m))
2393                   ENDIF
2394                ENDDO
2395
2396             ENDDO
2397          ENDIF
2398
2399          DO l = 0, 3
2400
2401             DO  m = 1, surf_lsm_v(l)%ns
2402!
2403!--             Spectral albedos for vegetation/pavement/water surfaces
2404                DO  ind_type = 0, 2
2405                   IF ( surf_lsm_v(l)%albedo_type(ind_type,m) /= 0 )  THEN
2406                      surf_lsm_v(l)%aldif(ind_type,m) =                        &
2407                            albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
2408                      surf_lsm_v(l)%asdif(ind_type,m) =                        &
2409                            albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2410                      surf_lsm_v(l)%aldir(ind_type,m) =                        &
2411                            albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
2412                      surf_lsm_v(l)%asdir(ind_type,m) =                        &
2413                            albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2414                      surf_lsm_v(l)%albedo(ind_type,m) =                       &
2415                            albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_type,m))
2416                   ENDIF
2417                ENDDO
2418             ENDDO
2419!
2420!--          For urban surface only if albedo has not been already initialized
2421!--          in the urban-surface model via the ASCII file.
2422             IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
2423                DO  m = 1, surf_usm_v(l)%ns
2424!
2425!--                Spectral albedos for wall/green/window surfaces
2426                   DO  ind_type = 0, 2
2427                      IF ( surf_usm_v(l)%albedo_type(ind_type,m) /= 0 )  THEN
2428                         surf_usm_v(l)%aldif(ind_type,m) =                     &
2429                            albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
2430                         surf_usm_v(l)%asdif(ind_type,m) =                     &
2431                            albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2432                         surf_usm_v(l)%aldir(ind_type,m) =                     &
2433                            albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
2434                         surf_usm_v(l)%asdir(ind_type,m) =                     &
2435                            albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2436                         surf_usm_v(l)%albedo(ind_type,m) =                    &
2437                            albedo_pars(2,surf_usm_v(l)%albedo_type(ind_type,m))
2438                      ENDIF
2439                   ENDDO
2440
2441                ENDDO
2442             ENDIF
2443          ENDDO
2444!
2445!--       Level 3 initialization at grid points where albedo type is zero.
2446!--       This case, spectral albedos are taken from file if available
2447          IF ( albedo_pars_f%from_file )  THEN
2448!
2449!--          Horizontal
2450             DO  m = 1, surf_lsm_h%ns
2451                i = surf_lsm_h%i(m)
2452                j = surf_lsm_h%j(m)
2453!
2454!--             Spectral albedos for vegetation/pavement/water surfaces
2455                DO  ind_type = 0, 2
2456                   IF ( surf_lsm_h%albedo_type(ind_type,m) == 0 )  THEN
2457                      IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )&
2458                         surf_lsm_h%albedo(ind_type,m) =                       &
2459                                                albedo_pars_f%pars_xy(0,j,i)
2460                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2461                         surf_lsm_h%aldir(ind_type,m) =                        &
2462                                                albedo_pars_f%pars_xy(1,j,i)
2463                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2464                         surf_lsm_h%aldif(ind_type,m) =                        &
2465                                                albedo_pars_f%pars_xy(1,j,i)
2466                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
2467                         surf_lsm_h%asdir(ind_type,m) =                        &
2468                                                albedo_pars_f%pars_xy(2,j,i)
2469                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
2470                         surf_lsm_h%asdif(ind_type,m) =                        &
2471                                                albedo_pars_f%pars_xy(2,j,i)
2472                   ENDIF
2473                ENDDO
2474             ENDDO
2475!
2476!--          For urban surface only if albedo has not been already initialized
2477!--          in the urban-surface model via the ASCII file.
2478             IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
2479                DO  m = 1, surf_usm_h%ns
2480                   i = surf_usm_h%i(m)
2481                   j = surf_usm_h%j(m)
2482!
2483!--                Broadband albedos for wall/green/window surfaces
2484                   DO  ind_type = 0, 2
2485                      IF ( surf_usm_h%albedo_type(ind_type,m) == 0 )  THEN
2486                         IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )&
2487                            surf_usm_h%albedo(ind_type,m) =                       &
2488                                                albedo_pars_f%pars_xy(0,j,i)
2489                      ENDIF
2490                   ENDDO
2491!
2492!--                Spectral albedos especially for building wall surfaces
2493                   IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )  THEN
2494                      surf_usm_h%aldir(ind_veg_wall,m) =                       &
2495                                                albedo_pars_f%pars_xy(1,j,i)
2496                      surf_usm_h%aldif(ind_veg_wall,m) =                       &
2497                                                albedo_pars_f%pars_xy(1,j,i)
2498                   ENDIF
2499                   IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )  THEN
2500                      surf_usm_h%asdir(ind_veg_wall,m) =                       &
2501                                                albedo_pars_f%pars_xy(2,j,i)
2502                      surf_usm_h%asdif(ind_veg_wall,m) =                       &
2503                                                albedo_pars_f%pars_xy(2,j,i)
2504                   ENDIF
2505!
2506!--                Spectral albedos especially for building green surfaces
2507                   IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )  THEN
2508                      surf_usm_h%aldir(ind_pav_green,m) =                      &
2509                                                albedo_pars_f%pars_xy(3,j,i)
2510                      surf_usm_h%aldif(ind_pav_green,m) =                      &
2511                                                albedo_pars_f%pars_xy(3,j,i)
2512                   ENDIF
2513                   IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )  THEN
2514                      surf_usm_h%asdir(ind_pav_green,m) =                      &
2515                                                albedo_pars_f%pars_xy(4,j,i)
2516                      surf_usm_h%asdif(ind_pav_green,m) =                      &
2517                                                albedo_pars_f%pars_xy(4,j,i)
2518                   ENDIF
2519!
2520!--                Spectral albedos especially for building window surfaces
2521                   IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill )  THEN
2522                      surf_usm_h%aldir(ind_wat_win,m) =                        &
2523                                                albedo_pars_f%pars_xy(5,j,i)
2524                      surf_usm_h%aldif(ind_wat_win,m) =                        &
2525                                                albedo_pars_f%pars_xy(5,j,i)
2526                   ENDIF
2527                   IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill )  THEN
2528                      surf_usm_h%asdir(ind_wat_win,m) =                        &
2529                                                albedo_pars_f%pars_xy(6,j,i)
2530                      surf_usm_h%asdif(ind_wat_win,m) =                        &
2531                                                albedo_pars_f%pars_xy(6,j,i)
2532                   ENDIF
2533
2534                ENDDO
2535             ENDIF
2536!
2537!--          Vertical
2538             DO  l = 0, 3
2539                ioff = surf_lsm_v(l)%ioff
2540                joff = surf_lsm_v(l)%joff
2541
2542                DO  m = 1, surf_lsm_v(l)%ns
2543                   i = surf_lsm_v(l)%i(m)
2544                   j = surf_lsm_v(l)%j(m)
2545!
2546!--                Spectral albedos for vegetation/pavement/water surfaces
2547                   DO  ind_type = 0, 2
2548                      IF ( surf_lsm_v(l)%albedo_type(ind_type,m) == 0 )  THEN
2549                         IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /=        &
2550                              albedo_pars_f%fill )                             &
2551                            surf_lsm_v(l)%albedo(ind_type,m) =                 &
2552                                          albedo_pars_f%pars_xy(0,j+joff,i+ioff)
2553                         IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
2554                              albedo_pars_f%fill )                             &
2555                            surf_lsm_v(l)%aldir(ind_type,m) =                  &
2556                                          albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2557                         IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
2558                              albedo_pars_f%fill )                             &
2559                            surf_lsm_v(l)%aldif(ind_type,m) =                  &
2560                                          albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2561                         IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /=        &
2562                              albedo_pars_f%fill )                             &
2563                            surf_lsm_v(l)%asdir(ind_type,m) =                  &
2564                                          albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2565                         IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /=        &
2566                              albedo_pars_f%fill )                             &
2567                            surf_lsm_v(l)%asdif(ind_type,m) =                  &
2568                                          albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2569                      ENDIF
2570                   ENDDO
2571                ENDDO
2572!
2573!--             For urban surface only if albedo has not been already initialized
2574!--             in the urban-surface model via the ASCII file.
2575                IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
2576                   ioff = surf_usm_v(l)%ioff
2577                   joff = surf_usm_v(l)%joff
2578
2579                   DO  m = 1, surf_usm_v(l)%ns
2580                      i = surf_usm_v(l)%i(m)
2581                      j = surf_usm_v(l)%j(m)
2582!
2583!--                   Broadband albedos for wall/green/window surfaces
2584                      DO  ind_type = 0, 2
2585                         IF ( surf_usm_v(l)%albedo_type(ind_type,m) == 0 )  THEN
2586                            IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /=     &
2587                                 albedo_pars_f%fill )                          &
2588                               surf_usm_v(l)%albedo(ind_type,m) =              &
2589                                             albedo_pars_f%pars_xy(0,j+joff,i+ioff)
2590                         ENDIF
2591                      ENDDO
2592!
2593!--                   Spectral albedos especially for building wall surfaces
2594                      IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=           &
2595                           albedo_pars_f%fill )  THEN
2596                         surf_usm_v(l)%aldir(ind_veg_wall,m) =                 &
2597                                         albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2598                         surf_usm_v(l)%aldif(ind_veg_wall,m) =                 &
2599                                         albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2600                      ENDIF
2601                      IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /=           &
2602                           albedo_pars_f%fill )  THEN
2603                         surf_usm_v(l)%asdir(ind_veg_wall,m) =                 &
2604                                         albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2605                         surf_usm_v(l)%asdif(ind_veg_wall,m) =                 &
2606                                         albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2607                      ENDIF
2608!                     
2609!--                   Spectral albedos especially for building green surfaces
2610                      IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /=           &
2611                           albedo_pars_f%fill )  THEN
2612                         surf_usm_v(l)%aldir(ind_pav_green,m) =                &
2613                                         albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2614                         surf_usm_v(l)%aldif(ind_pav_green,m) =                &
2615                                         albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2616                      ENDIF
2617                      IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /=           &
2618                           albedo_pars_f%fill )  THEN
2619                         surf_usm_v(l)%asdir(ind_pav_green,m) =                &
2620                                         albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2621                         surf_usm_v(l)%asdif(ind_pav_green,m) =                &
2622                                         albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2623                      ENDIF
2624!                     
2625!--                   Spectral albedos especially for building window surfaces
2626                      IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /=           &
2627                           albedo_pars_f%fill )  THEN
2628                         surf_usm_v(l)%aldir(ind_wat_win,m) =                  &
2629                                         albedo_pars_f%pars_xy(5,j+joff,i+ioff)
2630                         surf_usm_v(l)%aldif(ind_wat_win,m) =                  &
2631                                         albedo_pars_f%pars_xy(5,j+joff,i+ioff)
2632                      ENDIF
2633                      IF ( albedo_pars_f%pars_xy(6,j+joff,i+ioff) /=           &
2634                           albedo_pars_f%fill )  THEN
2635                         surf_usm_v(l)%asdir(ind_wat_win,m) =                  &
2636                                         albedo_pars_f%pars_xy(6,j+joff,i+ioff)
2637                         surf_usm_v(l)%asdif(ind_wat_win,m) =                  &
2638                                         albedo_pars_f%pars_xy(6,j+joff,i+ioff)
2639                      ENDIF
2640                   ENDDO
2641                ENDIF
2642             ENDDO
2643
2644          ENDIF
2645
2646!
2647!--       Calculate initial values of current (cosine of) the zenith angle and
2648!--       whether the sun is up
2649          CALL calc_zenith
2650!
2651!--       readjust date and time to its initial value
2652          CALL init_date_and_time
2653!
2654!--       Calculate initial surface albedo for different surfaces
2655          IF ( .NOT. constant_albedo )  THEN
2656#if defined( __netcdf )
2657!
2658!--          Horizontally aligned natural and urban surfaces
2659             CALL calc_albedo( surf_lsm_h )
2660             CALL calc_albedo( surf_usm_h )
2661!
2662!--          Vertically aligned natural and urban surfaces
2663             DO  l = 0, 3
2664                CALL calc_albedo( surf_lsm_v(l) )
2665                CALL calc_albedo( surf_usm_v(l) )
2666             ENDDO
2667#endif
2668          ELSE
2669!
2670!--          Initialize sun-inclination independent spectral albedos
2671!--          Horizontal surfaces
2672             IF ( surf_lsm_h%ns > 0 )  THEN
2673                surf_lsm_h%rrtm_aldir = surf_lsm_h%aldir
2674                surf_lsm_h%rrtm_asdir = surf_lsm_h%asdir
2675                surf_lsm_h%rrtm_aldif = surf_lsm_h%aldif
2676                surf_lsm_h%rrtm_asdif = surf_lsm_h%asdif
2677             ENDIF
2678             IF ( surf_usm_h%ns > 0 )  THEN
2679                surf_usm_h%rrtm_aldir = surf_usm_h%aldir
2680                surf_usm_h%rrtm_asdir = surf_usm_h%asdir
2681                surf_usm_h%rrtm_aldif = surf_usm_h%aldif
2682                surf_usm_h%rrtm_asdif = surf_usm_h%asdif
2683             ENDIF
2684!
2685!--          Vertical surfaces
2686             DO  l = 0, 3
2687                IF ( surf_lsm_v(l)%ns > 0 )  THEN
2688                   surf_lsm_v(l)%rrtm_aldir = surf_lsm_v(l)%aldir
2689                   surf_lsm_v(l)%rrtm_asdir = surf_lsm_v(l)%asdir
2690                   surf_lsm_v(l)%rrtm_aldif = surf_lsm_v(l)%aldif
2691                   surf_lsm_v(l)%rrtm_asdif = surf_lsm_v(l)%asdif
2692                ENDIF
2693                IF ( surf_usm_v(l)%ns > 0 )  THEN
2694                   surf_usm_v(l)%rrtm_aldir = surf_usm_v(l)%aldir
2695                   surf_usm_v(l)%rrtm_asdir = surf_usm_v(l)%asdir
2696                   surf_usm_v(l)%rrtm_aldif = surf_usm_v(l)%aldif
2697                   surf_usm_v(l)%rrtm_asdif = surf_usm_v(l)%asdif
2698                ENDIF
2699             ENDDO
2700
2701          ENDIF
2702
2703!
2704!--       Allocate 3d arrays of radiative fluxes and heating rates
2705          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
2706             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2707             rad_sw_in = 0.0_wp
2708          ENDIF
2709
2710          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
2711             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2712          ENDIF
2713
2714          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
2715             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2716             rad_sw_out = 0.0_wp
2717          ENDIF
2718
2719          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
2720             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2721          ENDIF
2722
2723          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
2724             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2725             rad_sw_hr = 0.0_wp
2726          ENDIF
2727
2728          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
2729             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2730             rad_sw_hr_av = 0.0_wp
2731          ENDIF
2732
2733          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
2734             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2735             rad_sw_cs_hr = 0.0_wp
2736          ENDIF
2737
2738          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
2739             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2740             rad_sw_cs_hr_av = 0.0_wp
2741          ENDIF
2742
2743          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
2744             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2745             rad_lw_in = 0.0_wp
2746          ENDIF
2747
2748          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
2749             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2750          ENDIF
2751
2752          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
2753             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2754            rad_lw_out = 0.0_wp
2755          ENDIF
2756
2757          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
2758             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2759          ENDIF
2760
2761          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
2762             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2763             rad_lw_hr = 0.0_wp
2764          ENDIF
2765
2766          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
2767             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2768             rad_lw_hr_av = 0.0_wp
2769          ENDIF
2770
2771          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
2772             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2773             rad_lw_cs_hr = 0.0_wp
2774          ENDIF
2775
2776          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
2777             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2778             rad_lw_cs_hr_av = 0.0_wp
2779          ENDIF
2780
2781          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2782          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2783          rad_sw_cs_in  = 0.0_wp
2784          rad_sw_cs_out = 0.0_wp
2785
2786          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2787          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2788          rad_lw_cs_in  = 0.0_wp
2789          rad_lw_cs_out = 0.0_wp
2790
2791!
2792!--       Allocate 1-element array for surface temperature
2793!--       (RRTMG anticipates an array as passed argument).
2794          ALLOCATE ( rrtm_tsfc(1) )
2795!
2796!--       Allocate surface emissivity.
2797!--       Values will be given directly before calling rrtm_lw.
2798          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
2799
2800!
2801!--       Initialize RRTMG, before check if files are existent
2802          INQUIRE( FILE='rrtmg_lw.nc', EXIST=lw_exists )
2803          IF ( .NOT. lw_exists )  THEN
2804             message_string = 'Input file rrtmg_lw.nc' //                &
2805                            '&for rrtmg missing. ' // &
2806                            '&Please provide <jobname>_lsw file in the INPUT directory.'
2807             CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 )
2808          ENDIF         
2809          INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists )
2810          IF ( .NOT. sw_exists )  THEN
2811             message_string = 'Input file rrtmg_sw.nc' //                &
2812                            '&for rrtmg missing. ' // &
2813                            '&Please provide <jobname>_rsw file in the INPUT directory.'
2814             CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 )
2815          ENDIF         
2816         
2817          IF ( lw_radiation )  CALL rrtmg_lw_ini ( c_p )
2818          IF ( sw_radiation )  CALL rrtmg_sw_ini ( c_p )
2819         
2820!
2821!--       Set input files for RRTMG
2822          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
2823          IF ( .NOT. snd_exists )  THEN
2824             rrtm_input_file = "rrtmg_lw.nc"
2825          ENDIF
2826
2827!
2828!--       Read vertical layers for RRTMG from sounding data
2829!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
2830!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
2831!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
2832          CALL read_sounding_data
2833
2834!
2835!--       Read trace gas profiles from file. This routine provides
2836!--       the rrtm_ arrays (1:nzt_rad+1)
2837          CALL read_trace_gas_data
2838#endif
2839       ENDIF
2840
2841!
2842!--    Perform user actions if required
2843       CALL user_init_radiation
2844
2845!
2846!--    Calculate radiative fluxes at model start
2847       SELECT CASE ( TRIM( radiation_scheme ) )
2848
2849          CASE ( 'rrtmg' )
2850             CALL radiation_rrtmg
2851
2852          CASE ( 'clear-sky' )
2853             CALL radiation_clearsky
2854
2855          CASE ( 'constant' )
2856             CALL radiation_constant
2857
2858          CASE DEFAULT
2859
2860       END SELECT
2861
2862! readjust date and time to its initial value
2863       CALL init_date_and_time
2864
2865!
2866!--    Find all discretized apparent solar positions for radiation interaction.
2867       IF ( radiation_interactions )  CALL radiation_presimulate_solar_pos
2868
2869!
2870!--    If required, read or calculate and write out the SVF
2871       IF ( radiation_interactions .AND. read_svf)  THEN
2872!
2873!--       Read sky-view factors and further required data from file
2874          CALL radiation_read_svf()
2875
2876       ELSEIF ( radiation_interactions .AND. .NOT. read_svf)  THEN
2877!
2878!--       calculate SFV and CSF
2879          CALL radiation_calc_svf()
2880       ENDIF
2881
2882       IF ( radiation_interactions .AND. write_svf)  THEN
2883!
2884!--       Write svf, csf svfsurf and csfsurf data to file
2885          CALL radiation_write_svf()
2886       ENDIF
2887
2888!
2889!--    Adjust radiative fluxes. In case of urban and land surfaces, also
2890!--    call an initial interaction.
2891       IF ( radiation_interactions )  THEN
2892          CALL radiation_interaction
2893       ENDIF
2894
2895       IF ( debug_output )  CALL debug_message( 'radiation_init', 'end' )
2896
2897       RETURN !todo: remove, I don't see what we need this for here
2898
2899    END SUBROUTINE radiation_init
2900
2901
2902!------------------------------------------------------------------------------!
2903! Description:
2904! ------------
2905!> A simple clear sky radiation model
2906!------------------------------------------------------------------------------!
2907    SUBROUTINE radiation_clearsky
2908
2909
2910       IMPLICIT NONE
2911
2912       INTEGER(iwp) ::  l         !< running index for surface orientation
2913       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
2914       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
2915       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
2916       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
2917
2918       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
2919
2920!
2921!--    Calculate current zenith angle
2922       CALL calc_zenith
2923
2924!
2925!--    Calculate sky transmissivity
2926       sky_trans = 0.6_wp + 0.2_wp * cos_zenith
2927
2928!
2929!--    Calculate value of the Exner function at model surface
2930!
2931!--    In case averaged radiation is used, calculate mean temperature and
2932!--    liquid water mixing ratio at the urban-layer top.
2933       IF ( average_radiation ) THEN
2934          pt1   = 0.0_wp
2935          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
2936
2937          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
2938          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
2939
2940#if defined( __parallel )     
2941          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2942          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2943          IF ( ierr /= 0 ) THEN
2944              WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
2945              FLUSH(9)
2946          ENDIF
2947
2948          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
2949              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2950              IF ( ierr /= 0 ) THEN
2951                  WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
2952                  FLUSH(9)
2953              ENDIF
2954          ENDIF
2955#else
2956          pt1 = pt1_l
2957          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
2958#endif
2959
2960          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1
2961!
2962!--       Finally, divide by number of grid points
2963          pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
2964       ENDIF
2965!
2966!--    Call clear-sky calculation for each surface orientation.
2967!--    First, horizontal surfaces
2968       surf => surf_lsm_h
2969       CALL radiation_clearsky_surf
2970       surf => surf_usm_h
2971       CALL radiation_clearsky_surf
2972!
2973!--    Vertical surfaces
2974       DO  l = 0, 3
2975          surf => surf_lsm_v(l)
2976          CALL radiation_clearsky_surf
2977          surf => surf_usm_v(l)
2978          CALL radiation_clearsky_surf
2979       ENDDO
2980
2981       CONTAINS
2982
2983          SUBROUTINE radiation_clearsky_surf
2984
2985             IMPLICIT NONE
2986
2987             INTEGER(iwp) ::  i         !< index x-direction
2988             INTEGER(iwp) ::  j         !< index y-direction
2989             INTEGER(iwp) ::  k         !< index z-direction
2990             INTEGER(iwp) ::  m         !< running index for surface elements
2991
2992             IF ( surf%ns < 1 )  RETURN
2993
2994!
2995!--          Calculate radiation fluxes and net radiation (rad_net) assuming
2996!--          homogeneous urban radiation conditions.
2997             IF ( average_radiation ) THEN       
2998
2999                k = nz_urban_t
3000
3001                surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
3002                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
3003               
3004                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
3005
3006                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
3007                                    + (1.0_wp - emissivity_urb) * surf%rad_lw_in
3008
3009                surf%rad_net = surf%rad_sw_in - surf%rad_sw_out                &
3010                             + surf%rad_lw_in - surf%rad_lw_out
3011
3012                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
3013                                           * (t_rad_urb)**3
3014
3015!
3016!--          Calculate radiation fluxes and net radiation (rad_net) for each surface
3017!--          element.
3018             ELSE
3019
3020                DO  m = 1, surf%ns
3021                   i = surf%i(m)
3022                   j = surf%j(m)
3023                   k = surf%k(m)
3024
3025                   surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith
3026
3027!
3028!--                Weighted average according to surface fraction.
3029!--                ATTENTION: when radiation interactions are switched on the
3030!--                calculated fluxes below are not actually used as they are
3031!--                overwritten in radiation_interaction.
3032                   surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
3033                                          surf%albedo(ind_veg_wall,m)          &
3034                                        + surf%frac(ind_pav_green,m) *         &
3035                                          surf%albedo(ind_pav_green,m)         &
3036                                        + surf%frac(ind_wat_win,m)   *         &
3037                                          surf%albedo(ind_wat_win,m) )         &
3038                                        * surf%rad_sw_in(m)
3039
3040                   surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
3041                                          surf%emissivity(ind_veg_wall,m)      &
3042                                        + surf%frac(ind_pav_green,m) *         &
3043                                          surf%emissivity(ind_pav_green,m)     &
3044                                        + surf%frac(ind_wat_win,m)   *         &
3045                                          surf%emissivity(ind_wat_win,m)       &
3046                                        )                                      &
3047                                        * sigma_sb                             &
3048                                        * ( surf%pt_surface(m) * exner(nzb) )**4
3049
3050                   surf%rad_lw_out_change_0(m) =                               &
3051                                      ( surf%frac(ind_veg_wall,m)  *           &
3052                                        surf%emissivity(ind_veg_wall,m)        &
3053                                      + surf%frac(ind_pav_green,m) *           &
3054                                        surf%emissivity(ind_pav_green,m)       &
3055                                      + surf%frac(ind_wat_win,m)   *           &
3056                                        surf%emissivity(ind_wat_win,m)         &
3057                                      ) * 4.0_wp * sigma_sb                    &
3058                                      * ( surf%pt_surface(m) * exner(nzb) )** 3
3059
3060
3061                   IF ( bulk_cloud_model  .OR.  cloud_droplets  )  THEN
3062                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
3063                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
3064                   ELSE
3065                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt(k,j,i) * exner(k))**4
3066                   ENDIF
3067
3068                   surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
3069                                   + surf%rad_lw_in(m) - surf%rad_lw_out(m)
3070
3071                ENDDO
3072
3073             ENDIF
3074
3075!
3076!--          Fill out values in radiation arrays
3077             DO  m = 1, surf%ns
3078                i = surf%i(m)
3079                j = surf%j(m)
3080                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3081                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3082                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3083                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3084             ENDDO
3085 
3086          END SUBROUTINE radiation_clearsky_surf
3087
3088    END SUBROUTINE radiation_clearsky
3089
3090
3091!------------------------------------------------------------------------------!
3092! Description:
3093! ------------
3094!> This scheme keeps the prescribed net radiation constant during the run
3095!------------------------------------------------------------------------------!
3096    SUBROUTINE radiation_constant
3097
3098
3099       IMPLICIT NONE
3100
3101       INTEGER(iwp) ::  l         !< running index for surface orientation
3102
3103       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
3104       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
3105       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
3106       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
3107
3108       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine
3109
3110!
3111!--    In case averaged radiation is used, calculate mean temperature and
3112!--    liquid water mixing ratio at the urban-layer top.
3113       IF ( average_radiation ) THEN   
3114          pt1   = 0.0_wp
3115          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
3116
3117          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
3118          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
3119
3120#if defined( __parallel )     
3121          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3122          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3123          IF ( ierr /= 0 ) THEN
3124              WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
3125              FLUSH(9)
3126          ENDIF
3127          IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
3128             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3129             IF ( ierr /= 0 ) THEN
3130                 WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
3131                 FLUSH(9)
3132             ENDIF
3133          ENDIF
3134#else
3135          pt1 = pt1_l
3136          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
3137#endif
3138          IF ( bulk_cloud_model  .OR.  cloud_droplets )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1
3139!
3140!--       Finally, divide by number of grid points
3141          pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
3142       ENDIF
3143
3144!
3145!--    First, horizontal surfaces
3146       surf => surf_lsm_h
3147       CALL radiation_constant_surf
3148       surf => surf_usm_h
3149       CALL radiation_constant_surf
3150!
3151!--    Vertical surfaces
3152       DO  l = 0, 3
3153          surf => surf_lsm_v(l)
3154          CALL radiation_constant_surf
3155          surf => surf_usm_v(l)
3156          CALL radiation_constant_surf
3157       ENDDO
3158
3159       CONTAINS
3160
3161          SUBROUTINE radiation_constant_surf
3162
3163             IMPLICIT NONE
3164
3165             INTEGER(iwp) ::  i         !< index x-direction
3166             INTEGER(iwp) ::  ioff      !< offset between surface element and adjacent grid point along x
3167             INTEGER(iwp) ::  j         !< index y-direction
3168             INTEGER(iwp) ::  joff      !< offset between surface element and adjacent grid point along y
3169             INTEGER(iwp) ::  k         !< index z-direction
3170             INTEGER(iwp) ::  koff      !< offset between surface element and adjacent grid point along z
3171             INTEGER(iwp) ::  m         !< running index for surface elements
3172
3173             IF ( surf%ns < 1 )  RETURN
3174
3175!--          Calculate homogenoeus urban radiation fluxes
3176             IF ( average_radiation ) THEN
3177
3178                surf%rad_net = net_radiation
3179
3180                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
3181
3182                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
3183                                    + ( 1.0_wp - emissivity_urb )             & ! shouldn't be this a bulk value -- emissivity_urb?
3184                                    * surf%rad_lw_in
3185
3186                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
3187                                           * t_rad_urb**3
3188
3189                surf%rad_sw_in = ( surf%rad_net - surf%rad_lw_in               &
3190                                     + surf%rad_lw_out )                       &
3191                                     / ( 1.0_wp - albedo_urb )
3192
3193                surf%rad_sw_out =  albedo_urb * surf%rad_sw_in
3194
3195!
3196!--          Calculate radiation fluxes for each surface element
3197             ELSE
3198!
3199!--             Determine index offset between surface element and adjacent
3200!--             atmospheric grid point
3201                ioff = surf%ioff
3202                joff = surf%joff
3203                koff = surf%koff
3204
3205!
3206!--             Prescribe net radiation and estimate the remaining radiative fluxes
3207                DO  m = 1, surf%ns
3208                   i = surf%i(m)
3209                   j = surf%j(m)
3210                   k = surf%k(m)
3211
3212                   surf%rad_net(m) = net_radiation
3213
3214                   IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
3215                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
3216                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
3217                   ELSE
3218                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb *                 &
3219                                             ( pt(k,j,i) * exner(k) )**4
3220                   ENDIF
3221
3222!
3223!--                Weighted average according to surface fraction.
3224                   surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
3225                                          surf%emissivity(ind_veg_wall,m)      &
3226                                        + surf%frac(ind_pav_green,m) *         &
3227                                          surf%emissivity(ind_pav_green,m)     &
3228                                        + surf%frac(ind_wat_win,m)   *         &
3229                                          surf%emissivity(ind_wat_win,m)       &
3230                                        )                                      &
3231                                      * sigma_sb                               &
3232                                      * ( surf%pt_surface(m) * exner(nzb) )**4
3233
3234                   surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m)   &
3235                                       + surf%rad_lw_out(m) )                  &
3236                                       / ( 1.0_wp -                            &
3237                                          ( surf%frac(ind_veg_wall,m)  *       &
3238                                            surf%albedo(ind_veg_wall,m)        &
3239                                         +  surf%frac(ind_pav_green,m) *       &
3240                                            surf%albedo(ind_pav_green,m)       &
3241                                         +  surf%frac(ind_wat_win,m)   *       &
3242                                            surf%albedo(ind_wat_win,m) )       &
3243                                         )
3244
3245                   surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
3246                                          surf%albedo(ind_veg_wall,m)          &
3247                                        + surf%frac(ind_pav_green,m) *         &
3248                                          surf%albedo(ind_pav_green,m)         &
3249                                        + surf%frac(ind_wat_win,m)   *         &
3250                                          surf%albedo(ind_wat_win,m) )         &
3251                                      * surf%rad_sw_in(m)
3252
3253                ENDDO
3254
3255             ENDIF
3256
3257!
3258!--          Fill out values in radiation arrays
3259             DO  m = 1, surf%ns
3260                i = surf%i(m)
3261                j = surf%j(m)
3262                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3263                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3264                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3265                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3266             ENDDO
3267
3268          END SUBROUTINE radiation_constant_surf
3269         
3270
3271    END SUBROUTINE radiation_constant
3272
3273!------------------------------------------------------------------------------!
3274! Description:
3275! ------------
3276!> Header output for radiation model
3277!------------------------------------------------------------------------------!
3278    SUBROUTINE radiation_header ( io )
3279
3280
3281       IMPLICIT NONE
3282 
3283       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
3284   
3285
3286       
3287!
3288!--    Write radiation model header
3289       WRITE( io, 3 )
3290
3291       IF ( radiation_scheme == "constant" )  THEN
3292          WRITE( io, 4 ) net_radiation
3293       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
3294          WRITE( io, 5 )
3295       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
3296          WRITE( io, 6 )
3297          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
3298          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
3299       ENDIF
3300
3301       IF ( albedo_type_f%from_file  .OR.  vegetation_type_f%from_file  .OR.   &
3302            pavement_type_f%from_file  .OR.  water_type_f%from_file  .OR.      &
3303            building_type_f%from_file )  THEN
3304             WRITE( io, 13 )
3305       ELSE 
3306          IF ( albedo_type == 0 )  THEN
3307             WRITE( io, 7 ) albedo
3308          ELSE
3309             WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
3310          ENDIF
3311       ENDIF
3312       IF ( constant_albedo )  THEN
3313          WRITE( io, 9 )
3314       ENDIF
3315       
3316       WRITE( io, 12 ) dt_radiation
3317 
3318
3319 3 FORMAT (//' Radiation model information:'/                                  &
3320              ' ----------------------------'/)
3321 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
3322           // 'W/m**2')
3323 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
3324                   ' default)')
3325 6 FORMAT ('    --> RRTMG scheme is used')
3326 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
3327 8 FORMAT (/'    Albedo is set for land surface type: ', A)
3328 9 FORMAT (/'    --> Albedo is fixed during the run')
332910 FORMAT (/'    --> Longwave radiation is disabled')
333011 FORMAT (/'    --> Shortwave radiation is disabled.')
333112 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
333213 FORMAT (/'    Albedo is set individually for each xy-location, according ', &
3333                 'to given surface type.')
3334
3335
3336    END SUBROUTINE radiation_header
3337   
3338
3339!------------------------------------------------------------------------------!
3340! Description:
3341! ------------
3342!> Parin for &radiation_parameters for radiation model
3343!------------------------------------------------------------------------------!
3344    SUBROUTINE radiation_parin
3345
3346
3347       IMPLICIT NONE
3348
3349       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
3350       
3351       NAMELIST /radiation_par/   albedo, albedo_lw_dif, albedo_lw_dir,         &
3352                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
3353                                  constant_albedo, dt_radiation, emissivity,    &
3354                                  lw_radiation, max_raytracing_dist,            &
3355                                  min_irrf_value, mrt_geom_human,               &
3356                                  mrt_include_sw, mrt_nlevels,                  &
3357                                  mrt_skip_roof, net_radiation, nrefsteps,      &
3358                                  plant_lw_interact, rad_angular_discretization,&
3359                                  radiation_interactions_on, radiation_scheme,  &
3360                                  raytrace_discrete_azims,                      &
3361                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
3362                                  skip_time_do_radiation, surface_reflections,  &
3363                                  svfnorm_report_thresh, sw_radiation,          &
3364                                  unscheduled_radiation_calls
3365
3366   
3367       NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir,    &
3368                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
3369                                  constant_albedo, dt_radiation, emissivity,    &
3370                                  lw_radiation, max_raytracing_dist,            &
3371                                  min_irrf_value, mrt_geom_human,               &
3372                                  mrt_include_sw, mrt_nlevels,                  &
3373                                  mrt_skip_roof, net_radiation, nrefsteps,      &
3374                                  plant_lw_interact, rad_angular_discretization,&
3375                                  radiation_interactions_on, radiation_scheme,  &
3376                                  raytrace_discrete_azims,                      &
3377                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
3378                                  skip_time_do_radiation, surface_reflections,  &
3379                                  svfnorm_report_thresh, sw_radiation,          &
3380                                  unscheduled_radiation_calls
3381   
3382       line = ' '
3383       
3384!
3385!--    Try to find radiation model namelist
3386       REWIND ( 11 )
3387       line = ' '
3388       DO WHILE ( INDEX( line, '&radiation_parameters' ) == 0 )
3389          READ ( 11, '(A)', END=12 )  line
3390       ENDDO
3391       BACKSPACE ( 11 )
3392
3393!
3394!--    Read user-defined namelist
3395       READ ( 11, radiation_parameters, ERR = 10 )
3396
3397!
3398!--    Set flag that indicates that the radiation model is switched on
3399       radiation = .TRUE.
3400
3401       GOTO 14
3402
3403 10    BACKSPACE( 11 )
3404       READ( 11 , '(A)') line
3405       CALL parin_fail_message( 'radiation_parameters', line )
3406!
3407!--    Try to find old namelist
3408 12    REWIND ( 11 )
3409       line = ' '
3410       DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
3411          READ ( 11, '(A)', END=14 )  line
3412       ENDDO
3413       BACKSPACE ( 11 )
3414
3415!
3416!--    Read user-defined namelist
3417       READ ( 11, radiation_par, ERR = 13, END = 14 )
3418
3419       message_string = 'namelist radiation_par is deprecated and will be ' // &
3420                     'removed in near future. Please use namelist ' //         &
3421                     'radiation_parameters instead'
3422       CALL message( 'radiation_parin', 'PA0487', 0, 1, 0, 6, 0 )
3423
3424!
3425!--    Set flag that indicates that the radiation model is switched on
3426       radiation = .TRUE.
3427
3428       IF ( .NOT.  radiation_interactions_on  .AND.  surface_reflections )  THEN
3429          message_string = 'surface_reflections is allowed only when '      // &
3430               'radiation_interactions_on is set to TRUE'
3431          CALL message( 'radiation_parin', 'PA0293',1, 2, 0, 6, 0 )
3432       ENDIF
3433
3434       GOTO 14
3435
3436 13    BACKSPACE( 11 )
3437       READ( 11 , '(A)') line
3438       CALL parin_fail_message( 'radiation_par', line )
3439
3440 14    CONTINUE
3441       
3442    END SUBROUTINE radiation_parin
3443
3444
3445!------------------------------------------------------------------------------!
3446! Description:
3447! ------------
3448!> Implementation of the RRTMG radiation_scheme
3449!------------------------------------------------------------------------------!
3450    SUBROUTINE radiation_rrtmg
3451
3452#if defined ( __rrtmg )
3453       USE indices,                                                            &
3454           ONLY:  nbgp
3455
3456       USE particle_attributes,                                                &
3457           ONLY:  grid_particles, number_of_particles, particles, prt_count
3458
3459       IMPLICIT NONE
3460
3461
3462       INTEGER(iwp) ::  i, j, k, l, m, n !< loop indices
3463       INTEGER(iwp) ::  k_topo     !< topography top index
3464
3465       REAL(wp)     ::  nc_rad, &    !< number concentration of cloud droplets
3466                        s_r2,   &    !< weighted sum over all droplets with r^2
3467                        s_r3         !< weighted sum over all droplets with r^3
3468
3469       REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av
3470       REAL(wp), DIMENSION(0:0)     :: zenith   !< to provide indexed array
3471!
3472!--    Just dummy arguments
3473       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rrtm_lw_taucld_dum,          &
3474                                                  rrtm_lw_tauaer_dum,          &
3475                                                  rrtm_sw_taucld_dum,          &
3476                                                  rrtm_sw_ssacld_dum,          &
3477                                                  rrtm_sw_asmcld_dum,          &
3478                                                  rrtm_sw_fsfcld_dum,          &
3479                                                  rrtm_sw_tauaer_dum,          &
3480                                                  rrtm_sw_ssaaer_dum,          &
3481                                                  rrtm_sw_asmaer_dum,          &
3482                                                  rrtm_sw_ecaer_dum
3483
3484!
3485!--    Calculate current (cosine of) zenith angle and whether the sun is up
3486       CALL calc_zenith     
3487       zenith(0) = cos_zenith
3488!
3489!--    Calculate surface albedo. In case average radiation is applied,
3490!--    this is not required.
3491#if defined( __netcdf )
3492       IF ( .NOT. constant_albedo )  THEN
3493!
3494!--       Horizontally aligned default, natural and urban surfaces
3495          CALL calc_albedo( surf_lsm_h    )
3496          CALL calc_albedo( surf_usm_h    )
3497!
3498!--       Vertically aligned default, natural and urban surfaces
3499          DO  l = 0, 3
3500             CALL calc_albedo( surf_lsm_v(l) )
3501             CALL calc_albedo( surf_usm_v(l) )
3502          ENDDO
3503       ENDIF
3504#endif
3505
3506!
3507!--    Prepare input data for RRTMG
3508
3509!
3510!--    In case of large scale forcing with surface data, calculate new pressure
3511!--    profile. nzt_rad might be modified by these calls and all required arrays
3512!--    will then be re-allocated
3513       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
3514          CALL read_sounding_data
3515          CALL read_trace_gas_data
3516       ENDIF
3517
3518
3519       IF ( average_radiation ) THEN
3520
3521          rrtm_asdir(1)  = albedo_urb
3522          rrtm_asdif(1)  = albedo_urb
3523          rrtm_aldir(1)  = albedo_urb
3524          rrtm_aldif(1)  = albedo_urb
3525
3526          rrtm_emis = emissivity_urb
3527!
3528!--       Calculate mean pt profile. Actually, only one height level is required.
3529          CALL calc_mean_profile( pt, 4 )
3530          pt_av = hom(:, 1, 4, 0)
3531         
3532          IF ( humidity )  THEN
3533             CALL calc_mean_profile( q, 41 )
3534             q_av  = hom(:, 1, 41, 0)
3535          ENDIF
3536!
3537!--       Prepare profiles of temperature and H2O volume mixing ratio
3538          rrtm_tlev(0,nzb+1) = t_rad_urb
3539
3540          IF ( bulk_cloud_model )  THEN
3541
3542             CALL calc_mean_profile( ql, 54 )
3543             ! average ql is now in hom(:, 1, 54, 0)
3544             ql_av = hom(:, 1, 54, 0)
3545             
3546             DO k = nzb+1, nzt+1
3547                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
3548                                 )**.286_wp + lv_d_cp * ql_av(k)
3549                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k))
3550             ENDDO
3551          ELSE
3552             DO k = nzb+1, nzt+1
3553                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
3554                                 )**.286_wp
3555             ENDDO
3556
3557             IF ( humidity )  THEN
3558                DO k = nzb+1, nzt+1
3559                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q_av(k)
3560                ENDDO
3561             ELSE
3562                rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
3563             ENDIF
3564          ENDIF
3565
3566!
3567!--       Avoid temperature/humidity jumps at the top of the PALM domain by
3568!--       linear interpolation from nzt+2 to nzt+7. Jumps are induced by
3569!--       discrepancies between the values in the  domain and those above that
3570!--       are prescribed in RRTMG
3571          DO k = nzt+2, nzt+7
3572             rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
3573                           + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
3574                           / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
3575                           * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3576
3577             rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
3578                           + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
3579                           / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
3580                           * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3581
3582          ENDDO
3583
3584!--       Linear interpolate to zw grid. Loop reaches one level further up
3585!--       due to the staggered grid in RRTMG
3586          DO k = nzb+2, nzt+8
3587             rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
3588                                rrtm_tlay(0,k-1))                           &
3589                                / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
3590                                * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
3591          ENDDO
3592
3593
3594!
3595!--       Calculate liquid water path and cloud fraction for each column.
3596!--       Note that LWP is required in g/m2 instead of kg/kg m.
3597          rrtm_cldfr  = 0.0_wp
3598          rrtm_reliq  = 0.0_wp
3599          rrtm_cliqwp = 0.0_wp
3600          rrtm_icld   = 0
3601
3602          IF ( bulk_cloud_model )  THEN
3603             DO k = nzb+1, nzt+1
3604                rrtm_cliqwp(0,k) =  ql_av(k) * 1000._wp *                   &
3605                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
3606                                    * 100._wp / g
3607
3608                IF ( rrtm_cliqwp(0,k) > 0._wp )  THEN
3609                   rrtm_cldfr(0,k) = 1._wp
3610                   IF ( rrtm_icld == 0 )  rrtm_icld = 1
3611
3612!
3613!--                Calculate cloud droplet effective radius
3614                   rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql_av(k)         &
3615                                     * rho_surface                          &
3616                                     / ( 4.0_wp * pi * nc_const * rho_l )   &
3617                                     )**0.33333333333333_wp                 &
3618                                     * EXP( LOG( sigma_gc )**2 )
3619!
3620!--                Limit effective radius
3621                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
3622                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
3623                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
3624                   ENDIF
3625                ENDIF
3626             ENDDO
3627          ENDIF
3628
3629!
3630!--       Set surface temperature
3631          rrtm_tsfc = t_rad_urb
3632         
3633          IF ( lw_radiation )  THEN       
3634         
3635             CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
3636             rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
3637             rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
3638             rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
3639             rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
3640             rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
3641             rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
3642             rrtm_reliq      , rrtm_lw_tauaer,                               &
3643             rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
3644             rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
3645             rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
3646
3647!
3648!--          Save fluxes
3649             DO k = nzb, nzt+1
3650                rad_lw_in(k,:,:)  = rrtm_lwdflx(0,k)
3651                rad_lw_out(k,:,:) = rrtm_lwuflx(0,k)
3652             ENDDO
3653             rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
3654!
3655!--          Save heating rates (convert from K/d to K/h).
3656!--          Further, even though an aggregated radiation is computed, map
3657!--          signle-column profiles on top of any topography, in order to
3658!--          obtain correct near surface radiation heating/cooling rates.
3659             DO  i = nxl, nxr
3660                DO  j = nys, nyn
3661                   k_topo = get_topography_top_index_ji( j, i, 's' )
3662                   DO k = k_topo+1, nzt+1
3663                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo)  * d_hours_day
3664                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo) * d_hours_day
3665                   ENDDO
3666                ENDDO
3667             ENDDO
3668
3669          ENDIF
3670
3671          IF ( sw_radiation .AND. sun_up )  THEN
3672             CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld     , rrtm_iaer      ,&
3673             rrtm_play      , rrtm_plev     , rrtm_tlay     , rrtm_tlev      ,&
3674             rrtm_tsfc      , rrtm_h2ovmr   , rrtm_o3vmr    , rrtm_co2vmr    ,&
3675             rrtm_ch4vmr    , rrtm_n2ovmr   , rrtm_o2vmr    , rrtm_asdir     ,&
3676             rrtm_asdif     , rrtm_aldir    , rrtm_aldif    , zenith         ,&
3677             0.0_wp         , day_of_year   , solar_constant, rrtm_inflgsw   ,&
3678             rrtm_iceflgsw  , rrtm_liqflgsw , rrtm_cldfr    , rrtm_sw_taucld ,&
3679             rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp    ,&
3680             rrtm_cliqwp    , rrtm_reice    , rrtm_reliq    , rrtm_sw_tauaer ,&
3681             rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx    ,&
3682             rrtm_swdflx    , rrtm_swhr     , rrtm_swuflxc  , rrtm_swdflxc   ,&
3683             rrtm_swhrc     , rrtm_dirdflux , rrtm_difdflux )
3684 
3685!
3686!--          Save fluxes:
3687!--          - whole domain
3688             DO k = nzb, nzt+1
3689                rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
3690                rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
3691             ENDDO
3692!--          - direct and diffuse SW at urban-surface-layer (required by RTM)
3693             rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb)
3694             rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb)
3695
3696!
3697!--          Save heating rates (convert from K/d to K/s)
3698             DO k = nzb+1, nzt+1
3699                rad_sw_hr(k,:,:)     = rrtm_swhr(0,k)  * d_hours_day
3700                rad_sw_cs_hr(k,:,:)  = rrtm_swhrc(0,k) * d_hours_day
3701             ENDDO
3702!
3703!--       Solar radiation is zero during night
3704          ELSE
3705             rad_sw_in  = 0.0_wp
3706             rad_sw_out = 0.0_wp
3707             rad_sw_in_dir(:,:) = 0.0_wp
3708             rad_sw_in_diff(:,:) = 0.0_wp
3709          ENDIF
3710!
3711!--    RRTMG is called for each (j,i) grid point separately, starting at the
3712!--    highest topography level. Here no RTM is used since average_radiation is false
3713       ELSE
3714!
3715!--       Loop over all grid points
3716          DO i = nxl, nxr
3717             DO j = nys, nyn
3718
3719!
3720!--             Prepare profiles of temperature and H2O volume mixing ratio
3721                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
3722                   rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * exner(nzb)
3723                ENDDO
3724                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3725                   rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * exner(nzb)
3726                ENDDO
3727
3728
3729                IF ( bulk_cloud_model )  THEN
3730                   DO k = nzb+1, nzt+1
3731                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                    &
3732                                        + lv_d_cp * ql(k,j,i)
3733                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
3734                   ENDDO
3735                ELSEIF ( cloud_droplets )  THEN
3736                   DO k = nzb+1, nzt+1
3737                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                     &
3738                                        + lv_d_cp * ql(k,j,i)
3739                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 
3740                   ENDDO
3741                ELSE
3742                   DO k = nzb+1, nzt+1
3743                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)
3744                   ENDDO
3745
3746                   IF ( humidity )  THEN
3747                      DO k = nzb+1, nzt+1
3748                         rrtm_h2ovmr(0,k) =  mol_mass_air_d_wv * q(k,j,i)
3749                      ENDDO   
3750                   ELSE
3751                      rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
3752                   ENDIF
3753                ENDIF
3754
3755!
3756!--             Avoid temperature/humidity jumps at the top of the LES domain by
3757!--             linear interpolation from nzt+2 to nzt+7
3758                DO k = nzt+2, nzt+7
3759                   rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                         &
3760                                 + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
3761                                 / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
3762                                 * ( rrtm_play(0,k)     - rrtm_play(0,nzt+1) )
3763
3764                   rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                     &
3765                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
3766                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
3767                              * ( rrtm_play(0,k)       - rrtm_play(0,nzt+1) )
3768
3769                ENDDO
3770
3771!--             Linear interpolate to zw grid
3772                DO k = nzb+2, nzt+8
3773                   rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -     &
3774                                      rrtm_tlay(0,k-1))                        &
3775                                      / ( rrtm_play(0,k) - rrtm_play(0,k-1) )  &
3776                                      * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
3777                ENDDO
3778
3779
3780!
3781!--             Calculate liquid water path and cloud fraction for each column.
3782!--             Note that LWP is required in g/m2 instead of kg/kg m.
3783                rrtm_cldfr  = 0.0_wp
3784                rrtm_reliq  = 0.0_wp
3785                rrtm_cliqwp = 0.0_wp
3786                rrtm_icld   = 0
3787
3788                IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
3789                   DO k = nzb+1, nzt+1
3790                      rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *              &
3791                                          (rrtm_plev(0,k) - rrtm_plev(0,k+1))  &
3792                                          * 100.0_wp / g
3793
3794                      IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
3795                         rrtm_cldfr(0,k) = 1.0_wp
3796                         IF ( rrtm_icld == 0 )  rrtm_icld = 1
3797
3798!
3799!--                      Calculate cloud droplet effective radius
3800                         IF ( bulk_cloud_model )  THEN
3801!
3802!--                         Calculete effective droplet radius. In case of using
3803!--                         cloud_scheme = 'morrison' and a non reasonable number
3804!--                         of cloud droplets the inital aerosol number 
3805!--                         concentration is considered.
3806                            IF ( microphysics_morrison )  THEN
3807                               IF ( nc(k,j,i) > 1.0E-20_wp )  THEN
3808                                  nc_rad = nc(k,j,i)
3809                               ELSE
3810                                  nc_rad = na_init
3811                               ENDIF
3812                            ELSE
3813                               nc_rad = nc_const
3814                            ENDIF 
3815
3816                            rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
3817                                              * rho_surface                       &
3818                                              / ( 4.0_wp * pi * nc_rad * rho_l )  &
3819                                              )**0.33333333333333_wp              &
3820                                              * EXP( LOG( sigma_gc )**2 )
3821
3822                         ELSEIF ( cloud_droplets )  THEN
3823                            number_of_particles = prt_count(k,j,i)
3824
3825                            IF (number_of_particles <= 0)  CYCLE
3826                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
3827                            s_r2 = 0.0_wp
3828                            s_r3 = 0.0_wp
3829
3830                            DO  n = 1, number_of_particles
3831                               IF ( particles(n)%particle_mask )  THEN
3832                                  s_r2 = s_r2 + particles(n)%radius**2 *       &
3833                                         particles(n)%weight_factor
3834                                  s_r3 = s_r3 + particles(n)%radius**3 *       &
3835                                         particles(n)%weight_factor
3836                               ENDIF
3837                            ENDDO
3838
3839                            IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
3840
3841                         ENDIF
3842
3843!
3844!--                      Limit effective radius
3845                         IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
3846                            rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
3847                            rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
3848                        ENDIF
3849                      ENDIF
3850                   ENDDO
3851                ENDIF
3852
3853!
3854!--             Write surface emissivity and surface temperature at current
3855!--             surface element on RRTMG-shaped array.
3856!--             Please note, as RRTMG is a single column model, surface attributes
3857!--             are only obtained from horizontally aligned surfaces (for
3858!--             simplicity). Taking surface attributes from horizontal and
3859!--             vertical walls would lead to multiple solutions. 
3860!--             Moreover, for natural- and urban-type surfaces, several surface
3861!--             classes can exist at a surface element next to each other.
3862!--             To obtain bulk parameters, apply a weighted average for these
3863!--             surfaces.
3864                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
3865                   rrtm_emis = surf_lsm_h%frac(ind_veg_wall,m)  *              &
3866                               surf_lsm_h%emissivity(ind_veg_wall,m)  +        &
3867                               surf_lsm_h%frac(ind_pav_green,m) *              &
3868                               surf_lsm_h%emissivity(ind_pav_green,m) +        & 
3869                               surf_lsm_h%frac(ind_wat_win,m)   *              &
3870                               surf_lsm_h%emissivity(ind_wat_win,m)
3871                   rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
3872                ENDDO             
3873                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3874                   rrtm_emis = surf_usm_h%frac(ind_veg_wall,m)  *              &
3875                               surf_usm_h%emissivity(ind_veg_wall,m)  +        &
3876                               surf_usm_h%frac(ind_pav_green,m) *              &
3877                               surf_usm_h%emissivity(ind_pav_green,m) +        & 
3878                               surf_usm_h%frac(ind_wat_win,m)   *              &
3879                               surf_usm_h%emissivity(ind_wat_win,m)
3880                   rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb)
3881                ENDDO
3882!
3883!--             Obtain topography top index (lower bound of RRTMG)
3884                k_topo = get_topography_top_index_ji( j, i, 's' )
3885
3886                IF ( lw_radiation )  THEN
3887!
3888!--                Due to technical reasons, copy optical depth to dummy arguments
3889!--                which are allocated on the exact size as the rrtmg_lw is called.
3890!--                As one dimesion is allocated with zero size, compiler complains
3891!--                that rank of the array does not match that of the
3892!--                assumed-shaped arguments in the RRTMG library. In order to
3893!--                avoid this, write to dummy arguments and give pass the entire
3894!--                dummy array. Seems to be the only existing work-around. 
3895                   ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) )
3896                   ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) )
3897
3898                   rrtm_lw_taucld_dum =                                        &
3899                               rrtm_lw_taucld(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1)
3900                   rrtm_lw_tauaer_dum =                                        &
3901                               rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)
3902
3903                   CALL rrtmg_lw( 1,                                           &                                       
3904                                  nzt_rad-k_topo,                              &
3905                                  rrtm_icld,                                   &
3906                                  rrtm_idrv,                                   &
3907                                  rrtm_play(:,k_topo+1:nzt_rad+1),             &
3908                                  rrtm_plev(:,k_topo+1:nzt_rad+2),             &
3909                                  rrtm_tlay(:,k_topo+1:nzt_rad+1),             &
3910                                  rrtm_tlev(:,k_topo+1:nzt_rad+2),             &
3911                                  rrtm_tsfc,                                   &
3912                                  rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),           &
3913                                  rrtm_o3vmr(:,k_topo+1:nzt_rad+1),            &
3914                                  rrtm_co2vmr(:,k_topo+1:nzt_rad+1),           &
3915                                  rrtm_ch4vmr(:,k_topo+1:nzt_rad+1),           &
3916                                  rrtm_n2ovmr(:,k_topo+1:nzt_rad+1),           &
3917                                  rrtm_o2vmr(:,k_topo+1:nzt_rad+1),            &
3918                                  rrtm_cfc11vmr(:,k_topo+1:nzt_rad+1),         &
3919                                  rrtm_cfc12vmr(:,k_topo+1:nzt_rad+1),         &
3920                                  rrtm_cfc22vmr(:,k_topo+1:nzt_rad+1),         &
3921                                  rrtm_ccl4vmr(:,k_topo+1:nzt_rad+1),          &
3922                                  rrtm_emis,                                   &
3923                                  rrtm_inflglw,                                &
3924                                  rrtm_iceflglw,                               &
3925                                  rrtm_liqflglw,                               &
3926                                  rrtm_cldfr(:,k_topo+1:nzt_rad+1),            &
3927                                  rrtm_lw_taucld_dum,                          &
3928                                  rrtm_cicewp(:,k_topo+1:nzt_rad+1),           &
3929                                  rrtm_cliqwp(:,k_topo+1:nzt_rad+1),           &
3930                                  rrtm_reice(:,k_topo+1:nzt_rad+1),            & 
3931                                  rrtm_reliq(:,k_topo+1:nzt_rad+1),            &
3932                                  rrtm_lw_tauaer_dum,                          &
3933                                  rrtm_lwuflx(:,k_topo:nzt_rad+1),             &
3934                                  rrtm_lwdflx(:,k_topo:nzt_rad+1),             &
3935                                  rrtm_lwhr(:,k_topo+1:nzt_rad+1),             &
3936                                  rrtm_lwuflxc(:,k_topo:nzt_rad+1),            &
3937                                  rrtm_lwdflxc(:,k_topo:nzt_rad+1),            &
3938                                  rrtm_lwhrc(:,k_topo+1:nzt_rad+1),            &
3939                                  rrtm_lwuflx_dt(:,k_topo:nzt_rad+1),          &
3940                                  rrtm_lwuflxc_dt(:,k_topo:nzt_rad+1) )
3941
3942                   DEALLOCATE ( rrtm_lw_taucld_dum )
3943                   DEALLOCATE ( rrtm_lw_tauaer_dum )
3944!
3945!--                Save fluxes
3946                   DO k = k_topo, nzt+1
3947                      rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
3948                      rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
3949                   ENDDO
3950
3951!
3952!--                Save heating rates (convert from K/d to K/h)
3953                   DO k = k_topo+1, nzt+1
3954                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo)  * d_hours_day
3955                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo) * d_hours_day
3956                   ENDDO
3957
3958!
3959!--                Save surface radiative fluxes and change in LW heating rate
3960!--                onto respective surface elements
3961!--                Horizontal surfaces
3962                   DO  m = surf_lsm_h%start_index(j,i),                        &
3963                           surf_lsm_h%end_index(j,i)
3964                      surf_lsm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
3965                      surf_lsm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
3966                      surf_lsm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
3967                   ENDDO             
3968                   DO  m = surf_usm_h%start_index(j,i),                        &
3969                           surf_usm_h%end_index(j,i)
3970                      surf_usm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
3971                      surf_usm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
3972                      surf_usm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
3973                   ENDDO
3974!
3975!--                Vertical surfaces. Fluxes are obtain at vertical level of the
3976!--                respective surface element
3977                   DO  l = 0, 3
3978                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
3979                              surf_lsm_v(l)%end_index(j,i)
3980                         k                                    = surf_lsm_v(l)%k(m)
3981                         surf_lsm_v(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k)
3982                         surf_lsm_v(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k)
3983                         surf_lsm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
3984                      ENDDO             
3985                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
3986                              surf_usm_v(l)%end_index(j,i)
3987                         k                                    = surf_usm_v(l)%k(m)
3988                         surf_usm_v(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k)
3989                         surf_usm_v(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k)
3990                         surf_usm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
3991                      ENDDO
3992                   ENDDO
3993
3994                ENDIF
3995
3996                IF ( sw_radiation .AND. sun_up )  THEN
3997!
3998!--                Get albedo for direct/diffusive long/shortwave radiation at
3999!--                current (y,x)-location from surface variables.
4000!--                Only obtain it from horizontal surfaces, as RRTMG is a single
4001!--                column model
4002!--                (Please note, only one loop will entered, controlled by
4003!--                start-end index.)
4004                   DO  m = surf_lsm_h%start_index(j,i),                        &
4005                           surf_lsm_h%end_index(j,i)
4006                      rrtm_asdir(1)  = SUM( surf_lsm_h%frac(:,m) *             &
4007                                            surf_lsm_h%rrtm_asdir(:,m) )
4008                      rrtm_asdif(1)  = SUM( surf_lsm_h%frac(:,m) *             &
4009                                            surf_lsm_h%rrtm_asdif(:,m) )
4010                      rrtm_aldir(1)  = SUM( surf_lsm_h%frac(:,m) *             &
4011                                            surf_lsm_h%rrtm_aldir(:,m) )
4012                      rrtm_aldif(1)  = SUM( surf_lsm_h%frac(:,m) *             &
4013                                            surf_lsm_h%rrtm_aldif(:,m) )
4014                   ENDDO             
4015                   DO  m = surf_usm_h%start_index(j,i),                        &
4016                           surf_usm_h%end_index(j,i)
4017                      rrtm_asdir(1)  = SUM( surf_usm_h%frac(:,m) *             &
4018                                            surf_usm_h%rrtm_asdir(:,m) )
4019                      rrtm_asdif(1)  = SUM( surf_usm_h%frac(:,m) *             &
4020                                            surf_usm_h%rrtm_asdif(:,m) )
4021                      rrtm_aldir(1)  = SUM( surf_usm_h%frac(:,m) *             &
4022                                            surf_usm_h%rrtm_aldir(:,m) )
4023                      rrtm_aldif(1)  = SUM( surf_usm_h%frac(:,m) *             &
4024                                            surf_usm_h%rrtm_aldif(:,m) )
4025                   ENDDO
4026!
4027!--                Due to technical reasons, copy optical depths and other
4028!--                to dummy arguments which are allocated on the exact size as the
4029!--                rrtmg_sw is called.
4030!--                As one dimesion is allocated with zero size, compiler complains
4031!--                that rank of the array does not match that of the
4032!--                assumed-shaped arguments in the RRTMG library. In order to
4033!--                avoid this, write to dummy arguments and give pass the entire
4034!--                dummy array. Seems to be the only existing work-around. 
4035                   ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
4036                   ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
4037                   ALLOCATE( rrtm_sw_asmcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
4038                   ALLOCATE( rrtm_sw_fsfcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
4039                   ALLOCATE( rrtm_sw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
4040                   ALLOCATE( rrtm_sw_ssaaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
4041                   ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
4042                   ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1)  )
4043     
4044                   rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
4045                   rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
4046                   rrtm_sw_asmcld_dum = rrtm_sw_asmcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
4047                   rrtm_sw_fsfcld_dum = rrtm_sw_fsfcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
4048                   rrtm_sw_tauaer_dum = rrtm_sw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
4049                   rrtm_sw_ssaaer_dum = rrtm_sw_ssaaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
4050                   rrtm_sw_asmaer_dum = rrtm_sw_asmaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
4051                   rrtm_sw_ecaer_dum  = rrtm_sw_ecaer(0:0,k_topo+1:nzt_rad+1,1:naerec+1)
4052
4053                   CALL rrtmg_sw( 1,                                           &
4054                                  nzt_rad-k_topo,                              &
4055                                  rrtm_icld,                                   &
4056                                  rrtm_iaer,                                   &
4057                                  rrtm_play(:,k_topo+1:nzt_rad+1),             &
4058                                  rrtm_plev(:,k_topo+1:nzt_rad+2),             &
4059                                  rrtm_tlay(:,k_topo+1:nzt_rad+1),             &
4060                                  rrtm_tlev(:,k_topo+1:nzt_rad+2),             &
4061                                  rrtm_tsfc,                                   &
4062                                  rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),           &                               
4063                                  rrtm_o3vmr(:,k_topo+1:nzt_rad+1),            &       
4064                                  rrtm_co2vmr(:,k_topo+1:nzt_rad+1),           &
4065                                  rrtm_ch4vmr(:,k_topo+1:nzt_rad+1),           &
4066                                  rrtm_n2ovmr(:,k_topo+1:nzt_rad+1),           &
4067                                  rrtm_o2vmr(:,k_topo+1:nzt_rad+1),            &
4068                                  rrtm_asdir,                                  & 
4069                                  rrtm_asdif,                                  &
4070                                  rrtm_aldir,                                  &
4071                                  rrtm_aldif,                                  &
4072                                  zenith,                                      &
4073                                  0.0_wp,                                      &
4074                                  day_of_year,                                 &
4075                                  solar_constant,                              &
4076                                  rrtm_inflgsw,                                &
4077                                  rrtm_iceflgsw,                               &
4078                                  rrtm_liqflgsw,                               &
4079                                  rrtm_cldfr(:,k_topo+1:nzt_rad+1),            &
4080                                  rrtm_sw_taucld_dum,                          &
4081                                  rrtm_sw_ssacld_dum,                          &
4082                                  rrtm_sw_asmcld_dum,                          &
4083                                  rrtm_sw_fsfcld_dum,                          &
4084                                  rrtm_cicewp(:,k_topo+1:nzt_rad+1),           &
4085                                  rrtm_cliqwp(:,k_topo+1:nzt_rad+1),           &
4086                                  rrtm_reice(:,k_topo+1:nzt_rad+1),            &
4087                                  rrtm_reliq(:,k_topo+1:nzt_rad+1),            &
4088                                  rrtm_sw_tauaer_dum,                          &
4089                                  rrtm_sw_ssaaer_dum,                          &
4090                                  rrtm_sw_asmaer_dum,                          &
4091                                  rrtm_sw_ecaer_dum,                           &
4092                                  rrtm_swuflx(:,k_topo:nzt_rad+1),             & 
4093                                  rrtm_swdflx(:,k_topo:nzt_rad+1),             & 
4094                                  rrtm_swhr(:,k_topo+1:nzt_rad+1),             & 
4095                                  rrtm_swuflxc(:,k_topo:nzt_rad+1),            & 
4096                                  rrtm_swdflxc(:,k_topo:nzt_rad+1),            &
4097                                  rrtm_swhrc(:,k_topo+1:nzt_rad+1),            &
4098                                  rrtm_dirdflux(:,k_topo:nzt_rad+1),           &
4099                                  rrtm_difdflux(:,k_topo:nzt_rad+1) )
4100
4101                   DEALLOCATE( rrtm_sw_taucld_dum )
4102                   DEALLOCATE( rrtm_sw_ssacld_dum )
4103                   DEALLOCATE( rrtm_sw_asmcld_dum )
4104                   DEALLOCATE( rrtm_sw_fsfcld_dum )
4105                   DEALLOCATE( rrtm_sw_tauaer_dum )
4106                   DEALLOCATE( rrtm_sw_ssaaer_dum )
4107                   DEALLOCATE( rrtm_sw_asmaer_dum )
4108                   DEALLOCATE( rrtm_sw_ecaer_dum )
4109!
4110!--                Save fluxes
4111                   DO k = nzb, nzt+1
4112                      rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
4113                      rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
4114                   ENDDO
4115!
4116!--                Save heating rates (convert from K/d to K/s)
4117                   DO k = nzb+1, nzt+1
4118                      rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
4119                      rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
4120                   ENDDO
4121
4122!
4123!--                Save surface radiative fluxes onto respective surface elements
4124!--                Horizontal surfaces
4125                   DO  m = surf_lsm_h%start_index(j,i),                        &
4126                           surf_lsm_h%end_index(j,i)
4127                      surf_lsm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
4128                      surf_lsm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
4129                   ENDDO             
4130                   DO  m = surf_usm_h%start_index(j,i),                        &
4131                           surf_usm_h%end_index(j,i)
4132                      surf_usm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
4133                      surf_usm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
4134                   ENDDO
4135!
4136!--                Vertical surfaces. Fluxes are obtain at respective vertical
4137!--                level of the surface element
4138                   DO  l = 0, 3
4139                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
4140                              surf_lsm_v(l)%end_index(j,i)
4141                         k                           = surf_lsm_v(l)%k(m)
4142                         surf_lsm_v(l)%rad_sw_in(m)  = rrtm_swdflx(0,k)
4143                         surf_lsm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k)
4144                      ENDDO             
4145                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
4146                              surf_usm_v(l)%end_index(j,i)
4147                         k                           = surf_usm_v(l)%k(m)
4148                         surf_usm_v(l)%rad_sw_in(m)  = rrtm_swdflx(0,k)
4149                         surf_usm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k)
4150                      ENDDO
4151                   ENDDO
4152!
4153!--             Solar radiation is zero during night
4154                ELSE
4155                   rad_sw_in  = 0.0_wp
4156                   rad_sw_out = 0.0_wp
4157!--             !!!!!!!! ATTENSION !!!!!!!!!!!!!!!
4158!--             Surface radiative fluxes should be also set to zero here                 
4159!--                Save surface radiative fluxes onto respective surface elements
4160!--                Horizontal surfaces
4161                   DO  m = surf_lsm_h%start_index(j,i),                        &
4162                           surf_lsm_h%end_index(j,i)
4163                      surf_lsm_h%rad_sw_in(m)     = 0.0_wp
4164                      surf_lsm_h%rad_sw_out(m)    = 0.0_wp
4165                   ENDDO             
4166                   DO  m = surf_usm_h%start_index(j,i),                        &
4167                           surf_usm_h%end_index(j,i)
4168                      surf_usm_h%rad_sw_in(m)     = 0.0_wp
4169                      surf_usm_h%rad_sw_out(m)    = 0.0_wp
4170                   ENDDO
4171!
4172!--                Vertical surfaces. Fluxes are obtain at respective vertical
4173!--                level of the surface element
4174                   DO  l = 0, 3
4175                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
4176                              surf_lsm_v(l)%end_index(j,i)
4177                         k                           = surf_lsm_v(l)%k(m)
4178                         surf_lsm_v(l)%rad_sw_in(m)  = 0.0_wp
4179                         surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp
4180                      ENDDO             
4181                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
4182                              surf_usm_v(l)%end_index(j,i)
4183                         k                           = surf_usm_v(l)%k(m)
4184                         surf_usm_v(l)%rad_sw_in(m)  = 0.0_wp
4185                         surf_usm_v(l)%rad_sw_out(m) = 0.0_wp
4186                      ENDDO
4187                   ENDDO
4188                ENDIF
4189
4190             ENDDO
4191          ENDDO
4192
4193       ENDIF
4194!
4195!--    Finally, calculate surface net radiation for surface elements.
4196       IF (  .NOT.  radiation_interactions  ) THEN
4197!--       First, for horizontal surfaces   
4198          DO  m = 1, surf_lsm_h%ns
4199             surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                   &
4200                                   - surf_lsm_h%rad_sw_out(m)                  &
4201                                   + surf_lsm_h%rad_lw_in(m)                   &
4202                                   - surf_lsm_h%rad_lw_out(m)
4203          ENDDO
4204          DO  m = 1, surf_usm_h%ns
4205             surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m)                   &
4206                                   - surf_usm_h%rad_sw_out(m)                  &
4207                                   + surf_usm_h%rad_lw_in(m)                   &
4208                                   - surf_usm_h%rad_lw_out(m)
4209          ENDDO
4210!
4211!--       Vertical surfaces.
4212!--       Todo: weight with azimuth and zenith angle according to their orientation!
4213          DO  l = 0, 3     
4214             DO  m = 1, surf_lsm_v(l)%ns
4215                surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)          &
4216                                         - surf_lsm_v(l)%rad_sw_out(m)         &
4217                                         + surf_lsm_v(l)%rad_lw_in(m)          &
4218                                         - surf_lsm_v(l)%rad_lw_out(m)
4219             ENDDO
4220             DO  m = 1, surf_usm_v(l)%ns
4221                surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)          &
4222                                         - surf_usm_v(l)%rad_sw_out(m)         &
4223                                         + surf_usm_v(l)%rad_lw_in(m)          &
4224                                         - surf_usm_v(l)%rad_lw_out(m)
4225             ENDDO
4226          ENDDO
4227       ENDIF
4228
4229
4230       CALL exchange_horiz( rad_lw_in,  nbgp )
4231       CALL exchange_horiz( rad_lw_out, nbgp )
4232       CALL exchange_horiz( rad_lw_hr,    nbgp )
4233       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
4234
4235       CALL exchange_horiz( rad_sw_in,  nbgp )
4236       CALL exchange_horiz( rad_sw_out, nbgp ) 
4237       CALL exchange_horiz( rad_sw_hr,    nbgp )
4238       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
4239
4240#endif
4241
4242    END SUBROUTINE radiation_rrtmg
4243
4244
4245!------------------------------------------------------------------------------!
4246! Description:
4247! ------------
4248!> Calculate the cosine of the zenith angle (variable is called zenith)
4249!------------------------------------------------------------------------------!
4250    SUBROUTINE calc_zenith
4251
4252       IMPLICIT NONE
4253
4254       REAL(wp) ::  declination,  & !< solar declination angle
4255                    hour_angle      !< solar hour angle
4256!
4257!--    Calculate current day and time based on the initial values and simulation
4258!--    time
4259       CALL calc_date_and_time
4260
4261!
4262!--    Calculate solar declination and hour angle   
4263       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day_of_year, KIND=wp) - decl_3) )
4264       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
4265
4266!
4267!--    Calculate cosine of solar zenith angle
4268       cos_zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
4269                                            * COS(hour_angle)
4270       cos_zenith = MAX(0.0_wp,cos_zenith)
4271
4272!
4273!--    Calculate solar directional vector
4274       IF ( sun_direction )  THEN
4275
4276!
4277!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
4278          sun_dir_lon = -SIN(hour_angle) * COS(declination)
4279
4280!
4281!--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
4282          sun_dir_lat = SIN(declination) * COS(lat) - COS(hour_angle) &
4283                              * COS(declination) * SIN(lat)
4284       ENDIF
4285
4286!
4287!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
4288       IF ( cos_zenith > 0.0_wp )  THEN
4289          sun_up = .TRUE.
4290       ELSE
4291          sun_up = .FALSE.
4292       END IF
4293
4294    END SUBROUTINE calc_zenith
4295
4296#if defined ( __rrtmg ) && defined ( __netcdf )
4297!------------------------------------------------------------------------------!
4298! Description:
4299! ------------
4300!> Calculates surface albedo components based on Briegleb (1992) and
4301!> Briegleb et al. (1986)
4302!------------------------------------------------------------------------------!
4303    SUBROUTINE calc_albedo( surf )
4304
4305        IMPLICIT NONE
4306
4307        INTEGER(iwp)    ::  ind_type !< running index surface tiles
4308        INTEGER(iwp)    ::  m        !< running index surface elements
4309
4310        TYPE(surf_type) ::  surf !< treated surfaces
4311
4312        IF ( sun_up  .AND.  .NOT. average_radiation )  THEN
4313
4314           DO  m = 1, surf%ns
4315!
4316!--           Loop over surface elements
4317              DO  ind_type = 0, SIZE( surf%albedo_type, 1 ) - 1
4318           
4319!
4320!--              Ocean
4321                 IF ( surf%albedo_type(ind_type,m) == 1 )  THEN
4322                    surf%rrtm_aldir(ind_type,m) = 0.026_wp /                    &
4323                                                ( cos_zenith**1.7_wp + 0.065_wp )&
4324                                     + 0.15_wp * ( cos_zenith - 0.1_wp )         &
4325                                               * ( cos_zenith - 0.5_wp )         &
4326                                               * ( cos_zenith - 1.0_wp )
4327                    surf%rrtm_asdir(ind_type,m) = surf%rrtm_aldir(ind_type,m)
4328!
4329!--              Snow
4330                 ELSEIF ( surf%albedo_type(ind_type,m) == 16 )  THEN
4331                    IF ( cos_zenith < 0.5_wp )  THEN
4332                       surf%rrtm_aldir(ind_type,m) =                           &
4333                                 0.5_wp * ( 1.0_wp - surf%aldif(ind_type,m) )  &
4334                                        * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
4335                                        * cos_zenith ) ) - 1.0_wp
4336                       surf%rrtm_asdir(ind_type,m) =                           &
4337                                 0.5_wp * ( 1.0_wp - surf%asdif(ind_type,m) )  &
4338                                        * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
4339                                        * cos_zenith ) ) - 1.0_wp
4340
4341                       surf%rrtm_aldir(ind_type,m) =                           &
4342                                       MIN(0.98_wp, surf%rrtm_aldir(ind_type,m))
4343                       surf%rrtm_asdir(ind_type,m) =                           &
4344                                       MIN(0.98_wp, surf%rrtm_asdir(ind_type,m))
4345                    ELSE
4346                       surf%rrtm_aldir(ind_type,m) = surf%aldif(ind_type,m)
4347                       surf%rrtm_asdir(ind_type,m) = surf%asdif(ind_type,m)
4348                    ENDIF
4349!
4350!--              Sea ice
4351                 ELSEIF ( surf%albedo_type(ind_type,m) == 15 )  THEN
4352                    surf%rrtm_aldir(ind_type,m) = surf%aldif(ind_type,m)
4353                    surf%rrtm_asdir(ind_type,m) = surf%asdif(ind_type,m)
4354
4355!
4356!--              Asphalt
4357                 ELSEIF ( surf%albedo_type(ind_type,m) == 17 )  THEN
4358                    surf%rrtm_aldir(ind_type,m) = surf%aldif(ind_type,m)
4359                    surf%rrtm_asdir(ind_type,m) = surf%asdif(ind_type,m)
4360
4361
4362!
4363!--              Bare soil
4364                 ELSEIF ( surf%albedo_type(ind_type,m) == 18 )  THEN
4365                    surf%rrtm_aldir(ind_type,m) = surf%aldif(ind_type,m)
4366                    surf%rrtm_asdir(ind_type,m) = surf%asdif(ind_type,m)
4367
4368!
4369!--              Land surfaces
4370                 ELSE
4371                    SELECT CASE ( surf%albedo_type(ind_type,m) )
4372
4373!
4374!--                    Surface types with strong zenith dependence
4375                       CASE ( 1, 2, 3, 4, 11, 12, 13 )
4376                          surf%rrtm_aldir(ind_type,m) =                        &
4377                                surf%aldif(ind_type,m) * 1.4_wp /              &
4378                                           ( 1.0_wp + 0.8_wp * cos_zenith )
4379                          surf%rrtm_asdir(ind_type,m) =                        &
4380                                surf%asdif(ind_type,m) * 1.4_wp /              &
4381                                           ( 1.0_wp + 0.8_wp * cos_zenith )
4382!
4383!--                    Surface types with weak zenith dependence
4384                       CASE ( 5, 6, 7, 8, 9, 10, 14 )
4385                          surf%rrtm_aldir(ind_type,m) =                        &
4386                                surf%aldif(ind_type,m) * 1.1_wp /              &
4387                                           ( 1.0_wp + 0.2_wp * cos_zenith )
4388                          surf%rrtm_asdir(ind_type,m) =                        &
4389                                surf%asdif(ind_type,m) * 1.1_wp /              &
4390                                           ( 1.0_wp + 0.2_wp * cos_zenith )
4391
4392                       CASE DEFAULT
4393
4394                    END SELECT
4395                 ENDIF
4396!
4397!--              Diffusive albedo is taken from Table 2
4398                 surf%rrtm_aldif(ind_type,m) = surf%aldif(ind_type,m)
4399                 surf%rrtm_asdif(ind_type,m) = surf%asdif(ind_type,m)
4400              ENDDO
4401           ENDDO
4402!
4403!--     Set albedo in case of average radiation
4404        ELSEIF ( sun_up  .AND.  average_radiation )  THEN
4405           surf%rrtm_asdir = albedo_urb
4406           surf%rrtm_asdif = albedo_urb
4407           surf%rrtm_aldir = albedo_urb
4408           surf%rrtm_aldif = albedo_urb 
4409!
4410!--     Darkness
4411        ELSE
4412           surf%rrtm_aldir = 0.0_wp
4413           surf%rrtm_asdir = 0.0_wp
4414           surf%rrtm_aldif = 0.0_wp
4415           surf%rrtm_asdif = 0.0_wp
4416        ENDIF
4417
4418    END SUBROUTINE calc_albedo
4419
4420!------------------------------------------------------------------------------!
4421! Description:
4422! ------------
4423!> Read sounding data (pressure and temperature) from RADIATION_DATA.
4424!------------------------------------------------------------------------------!
4425    SUBROUTINE read_sounding_data
4426
4427       IMPLICIT NONE
4428
4429       INTEGER(iwp) :: id,           & !< NetCDF id of input file
4430                       id_dim_zrad,  & !< pressure level id in the NetCDF file
4431                       id_var,       & !< NetCDF variable id
4432                       k,            & !< loop index
4433                       nz_snd,       & !< number of vertical levels in the sounding data
4434                       nz_snd_start, & !< start vertical index for sounding data to be used
4435                       nz_snd_end      !< end vertical index for souding data to be used
4436
4437       REAL(wp) :: t_surface           !< actual surface temperature
4438
4439       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
4440                                               t_snd_tmp      !< temporary temperature profile (sounding)
4441
4442!
4443!--    In case of updates, deallocate arrays first (sufficient to check one
4444!--    array as the others are automatically allocated). This is required
4445!--    because nzt_rad might change during the update
4446       IF ( ALLOCATED ( hyp_snd ) )  THEN
4447          DEALLOCATE( hyp_snd )
4448          DEALLOCATE( t_snd )
4449          DEALLOCATE ( rrtm_play )
4450          DEALLOCATE ( rrtm_plev )
4451          DEALLOCATE ( rrtm_tlay )
4452          DEALLOCATE ( rrtm_tlev )
4453
4454          DEALLOCATE ( rrtm_cicewp )
4455          DEALLOCATE ( rrtm_cldfr )
4456          DEALLOCATE ( rrtm_cliqwp )
4457          DEALLOCATE ( rrtm_reice )
4458          DEALLOCATE ( rrtm_reliq )
4459          DEALLOCATE ( rrtm_lw_taucld )
4460          DEALLOCATE ( rrtm_lw_tauaer )
4461
4462          DEALLOCATE ( rrtm_lwdflx  )
4463          DEALLOCATE ( rrtm_lwdflxc )
4464          DEALLOCATE ( rrtm_lwuflx  )
4465          DEALLOCATE ( rrtm_lwuflxc )
4466          DEALLOCATE ( rrtm_lwuflx_dt )
4467          DEALLOCATE ( rrtm_lwuflxc_dt )
4468          DEALLOCATE ( rrtm_lwhr  )
4469          DEALLOCATE ( rrtm_lwhrc )
4470
4471          DEALLOCATE ( rrtm_sw_taucld )
4472          DEALLOCATE ( rrtm_sw_ssacld )
4473          DEALLOCATE ( rrtm_sw_asmcld )
4474          DEALLOCATE ( rrtm_sw_fsfcld )
4475          DEALLOCATE ( rrtm_sw_tauaer )
4476          DEALLOCATE ( rrtm_sw_ssaaer )
4477          DEALLOCATE ( rrtm_sw_asmaer ) 
4478          DEALLOCATE ( rrtm_sw_ecaer )   
4479 
4480          DEALLOCATE ( rrtm_swdflx  )
4481          DEALLOCATE ( rrtm_swdflxc )
4482          DEALLOCATE ( rrtm_swuflx  )
4483          DEALLOCATE ( rrtm_swuflxc )
4484          DEALLOCATE ( rrtm_swhr  )
4485          DEALLOCATE ( rrtm_swhrc )
4486          DEALLOCATE ( rrtm_dirdflux )
4487          DEALLOCATE ( rrtm_difdflux )
4488
4489       ENDIF
4490
4491!
4492!--    Open file for reading
4493       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
4494       CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
4495
4496!
4497!--    Inquire dimension of z axis and save in nz_snd
4498       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
4499       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
4500       CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
4501
4502!
4503! !--    Allocate temporary array for storing pressure data
4504       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
4505       hyp_snd_tmp = 0.0_wp
4506
4507
4508!--    Read pressure from file
4509       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
4510       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
4511                               count = (/nz_snd/) )
4512       CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
4513
4514!
4515!--    Allocate temporary array for storing temperature data
4516       ALLOCATE( t_snd_tmp(1:nz_snd) )
4517       t_snd_tmp = 0.0_wp
4518
4519!
4520!--    Read temperature from file
4521       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
4522       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
4523                               count = (/nz_snd/) )
4524       CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
4525
4526!
4527!--    Calculate start of sounding data
4528       nz_snd_start = nz_snd + 1
4529       nz_snd_end   = nz_snd + 1
4530
4531!
4532!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
4533!--    in Pa, hyp_snd in hPa).
4534       DO  k = 1, nz_snd
4535          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
4536             nz_snd_start = k
4537             EXIT
4538          END IF
4539       END DO
4540
4541       IF ( nz_snd_start <= nz_snd )  THEN
4542          nz_snd_end = nz_snd
4543       END IF
4544
4545
4546!
4547!--    Calculate of total grid points for RRTMG calculations
4548       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
4549
4550!
4551!--    Save data above LES domain in hyp_snd, t_snd
4552       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
4553       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
4554       hyp_snd = 0.0_wp
4555       t_snd = 0.0_wp
4556
4557       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
4558       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start+1:nz_snd_end)
4559
4560       nc_stat = NF90_CLOSE( id )
4561
4562!
4563!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
4564!--    top of the LES domain. This routine does not consider horizontal or
4565!--    vertical variability of pressure and temperature
4566       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
4567       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
4568
4569       t_surface = pt_surface * exner(nzb)
4570       DO k = nzb+1, nzt+1
4571          rrtm_play(0,k) = hyp(k) * 0.01_wp
4572          rrtm_plev(0,k) = barometric_formula(zw(k-1),                         &
4573                              pt_surface * exner(nzb), &
4574                              surface_pressure )
4575       ENDDO
4576
4577       DO k = nzt+2, nzt_rad
4578          rrtm_play(0,k) = hyp_snd(k)
4579          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
4580       ENDDO
4581       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
4582                                   1.5 * hyp_snd(nzt_rad)                      &
4583                                 - 0.5 * hyp_snd(nzt_rad-1) )
4584       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
4585                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
4586
4587       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
4588
4589!
4590!--    Calculate temperature/humidity levels at top of the LES domain.
4591!--    Currently, the temperature is taken from sounding data (might lead to a
4592!--    temperature jump at interface. To do: Humidity is currently not
4593!--    calculated above the LES domain.
4594       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
4595       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
4596
4597       DO k = nzt+8, nzt_rad
4598          rrtm_tlay(0,k)   = t_snd(k)
4599       ENDDO
4600       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
4601                                - rrtm_tlay(0,nzt_rad-1)
4602       DO k = nzt+9, nzt_rad+1
4603          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
4604                             - rrtm_tlay(0,k-1))                               &
4605                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
4606                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
4607       ENDDO
4608
4609       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
4610                                  - rrtm_tlev(0,nzt_rad)
4611!
4612!--    Allocate remaining RRTMG arrays
4613       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
4614       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
4615       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
4616       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
4617       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
4618       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
4619       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
4620       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
4621       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
4622       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
4623       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
4624       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
4625       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
4626       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
4627       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
4628
4629!
4630!--    The ice phase is currently not considered in PALM
4631       rrtm_cicewp = 0.0_wp
4632       rrtm_reice  = 0.0_wp
4633
4634!
4635!--    Set other parameters (move to NAMELIST parameters in the future)
4636       rrtm_lw_tauaer = 0.0_wp
4637       rrtm_lw_taucld = 0.0_wp
4638       rrtm_sw_taucld = 0.0_wp
4639       rrtm_sw_ssacld = 0.0_wp
4640       rrtm_sw_asmcld = 0.0_wp
4641       rrtm_sw_fsfcld = 0.0_wp
4642       rrtm_sw_tauaer = 0.0_wp
4643       rrtm_sw_ssaaer = 0.0_wp
4644       rrtm_sw_asmaer = 0.0_wp
4645       rrtm_sw_ecaer  = 0.0_wp
4646
4647
4648       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
4649       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
4650       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
4651       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
4652       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
4653       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
4654       ALLOCATE ( rrtm_dirdflux(0:0,nzb:nzt_rad+1) )
4655       ALLOCATE ( rrtm_difdflux(0:0,nzb:nzt_rad+1) )
4656
4657       rrtm_swdflx  = 0.0_wp
4658       rrtm_swuflx  = 0.0_wp
4659       rrtm_swhr    = 0.0_wp 
4660       rrtm_swuflxc = 0.0_wp
4661       rrtm_swdflxc = 0.0_wp
4662       rrtm_swhrc   = 0.0_wp
4663       rrtm_dirdflux = 0.0_wp
4664       rrtm_difdflux = 0.0_wp
4665
4666       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
4667       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
4668       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
4669       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
4670       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
4671       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
4672
4673       rrtm_lwdflx  = 0.0_wp
4674       rrtm_lwuflx  = 0.0_wp
4675       rrtm_lwhr    = 0.0_wp 
4676       rrtm_lwuflxc = 0.0_wp
4677       rrtm_lwdflxc = 0.0_wp
4678       rrtm_lwhrc   = 0.0_wp
4679
4680       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
4681       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
4682
4683       rrtm_lwuflx_dt = 0.0_wp
4684       rrtm_lwuflxc_dt = 0.0_wp
4685
4686    END SUBROUTINE read_sounding_data
4687
4688
4689!------------------------------------------------------------------------------!
4690! Description:
4691! ------------
4692!> Read trace gas data from file and convert into trace gas paths / volume
4693!> mixing ratios. If a user-defined input file is provided it needs to follow
4694!> the convections used in RRTMG (see respective netCDF files shipped with
4695!> RRTMG)
4696!------------------------------------------------------------------------------!
4697    SUBROUTINE read_trace_gas_data
4698
4699       USE rrsw_ncpar
4700
4701       IMPLICIT NONE
4702
4703       INTEGER(iwp), PARAMETER :: num_trace_gases = 10 !< number of trace gases (absorbers)
4704
4705       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
4706           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
4707                           'CFC11', 'CFC12', 'CFC22', 'CCL4 ', 'H2O  '/)
4708
4709       INTEGER(iwp) :: id,     & !< NetCDF id
4710                       k,      & !< loop index
4711                       m,      & !< loop index
4712                       n,      & !< loop index
4713                       nabs,   & !< number of absorbers
4714                       np,     & !< number of pressure levels
4715                       id_abs, & !< NetCDF id of the respective absorber
4716                       id_dim, & !< NetCDF id of asborber's dimension
4717                       id_var    !< NetCDf id ot the absorber
4718
4719       REAL(wp) :: p_mls_l, &    !< pressure lower limit for interpolation
4720                   p_mls_u, &    !< pressure upper limit for interpolation
4721                   p_wgt_l, &    !< pressure weight lower limit for interpolation
4722                   p_wgt_u, &    !< pressure weight upper limit for interpolation
4723                   p_mls_m       !< mean pressure between upper and lower limits
4724
4725
4726       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,          & !< pressure levels for the absorbers
4727                                                 rrtm_play_tmp,  & !< temporary array for pressure zu-levels
4728                                                 rrtm_plev_tmp,  & !< temporary array for pressure zw-levels
4729                                                 trace_path_tmp    !< temporary array for storing trace gas path data
4730
4731       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
4732                                                 trace_mls_path, & !< array for storing trace gas path data
4733                                                 trace_mls_tmp     !< temporary array for storing trace gas data
4734
4735
4736!
4737!--    In case of updates, deallocate arrays first (sufficient to check one
4738!--    array as the others are automatically allocated)
4739       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
4740          DEALLOCATE ( rrtm_o3vmr  )
4741          DEALLOCATE ( rrtm_co2vmr )
4742          DEALLOCATE ( rrtm_ch4vmr )
4743          DEALLOCATE ( rrtm_n2ovmr )
4744          DEALLOCATE ( rrtm_o2vmr  )
4745          DEALLOCATE ( rrtm_cfc11vmr )
4746          DEALLOCATE ( rrtm_cfc12vmr )
4747          DEALLOCATE ( rrtm_cfc22vmr )
4748          DEALLOCATE ( rrtm_ccl4vmr  )
4749          DEALLOCATE ( rrtm_h2ovmr  )     
4750       ENDIF
4751
4752!
4753!--    Allocate trace gas profiles
4754       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
4755       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
4756       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
4757       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
4758       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
4759       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
4760       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
4761       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
4762       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
4763       ALLOCATE ( rrtm_h2ovmr(0:0,1:nzt_rad+1)  )
4764
4765!
4766!--    Open file for reading
4767       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
4768       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
4769!
4770!--    Inquire dimension ids and dimensions
4771       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
4772       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4773       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
4774       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4775
4776       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
4777       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4778       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
4779       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4780   
4781
4782!
4783!--    Allocate pressure, and trace gas arrays     
4784       ALLOCATE( p_mls(1:np) )
4785       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
4786       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
4787
4788
4789       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
4790       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4791       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
4792       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4793
4794       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
4795       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4796       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
4797       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
4798
4799
4800!
4801!--    Write absorber amounts (m