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

Last change on this file since 3943 was 3943, checked in by maronga, 5 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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