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

Last change on this file since 2550 was 2547, checked in by schwenkel, 7 years ago

extended by cloud_droplets option

  • Property svn:keywords set to Id
File size: 148.7 KB
Line 
1!> @file radiation_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
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 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: radiation_model_mod.f90 2547 2017-10-16 12:41:56Z boeske $
27! extended by cloud_droplets option, minor bugfix and correct calculation of
28! cloud droplet number concentration
29!
30! 2544 2017-10-13 18:09:32Z maronga
31! Moved date and time quantitis to separate module date_and_time_mod
32!
33! 2512 2017-10-04 08:26:59Z raasch
34! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
35! no output of ghost layer data
36!
37! 2504 2017-09-27 10:36:13Z maronga
38! Updates pavement types and albedo parameters
39!
40! 2328 2017-08-03 12:34:22Z maronga
41! Emissivity can now be set individually for each pixel.
42! Albedo type can be inferred from land surface model.
43! Added default albedo type for bare soil
44!
45! 2318 2017-07-20 17:27:44Z suehring
46! Get topography top index via Function call
47!
48! 2317 2017-07-20 17:27:19Z suehring
49! Improved syntax layout
50!
51! 2298 2017-06-29 09:28:18Z raasch
52! type of write_binary changed from CHARACTER to LOGICAL
53!
54! 2296 2017-06-28 07:53:56Z maronga
55! Added output of rad_sw_out for radiation_scheme = 'constant'
56!
57! 2270 2017-06-09 12:18:47Z maronga
58! Numbering changed (2 timeseries removed)
59!
60! 2249 2017-06-06 13:58:01Z sward
61! Allow for RRTMG runs without humidity/cloud physics
62!
63! 2248 2017-06-06 13:52:54Z sward
64! Error no changed
65!
66! 2233 2017-05-30 18:08:54Z suehring
67!
68! 2232 2017-05-30 17:47:52Z suehring
69! Adjustments to new topography concept
70! Bugfix in read restart
71!
72! 2200 2017-04-11 11:37:51Z suehring
73! Bugfix in call of exchange_horiz_2d and read restart data
74!
75! 2163 2017-03-01 13:23:15Z schwenkel
76! Bugfix in radiation_check_data_output
77!
78! 2157 2017-02-22 15:10:35Z suehring
79! Bugfix in read_restart data
80!
81! 2011 2016-09-19 17:29:57Z kanani
82! Removed CALL of auxiliary SUBROUTINE get_usm_info,
83! flag urban_surface is now defined in module control_parameters.
84!
85! 2007 2016-08-24 15:47:17Z kanani
86! Added calculation of solar directional vector for new urban surface
87! model,
88! accounted for urban_surface model in radiation_check_parameters,
89! correction of comments for zenith angle.
90!
91! 2000 2016-08-20 18:09:15Z knoop
92! Forced header and separation lines into 80 columns
93!
94! 1976 2016-07-27 13:28:04Z maronga
95! Output of 2D/3D/masked data is now directly done within this module. The
96! radiation schemes have been simplified for better usability so that
97! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
98! the radiation code used.
99!
100! 1856 2016-04-13 12:56:17Z maronga
101! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
102!
103! 1853 2016-04-11 09:00:35Z maronga
104! Added routine for radiation_scheme = constant.
105
106! 1849 2016-04-08 11:33:18Z hoffmann
107! Adapted for modularization of microphysics
108!
109! 1826 2016-04-07 12:01:39Z maronga
110! Further modularization.
111!
112! 1788 2016-03-10 11:01:04Z maronga
113! Added new albedo class for pavements / roads.
114!
115! 1783 2016-03-06 18:36:17Z raasch
116! palm-netcdf-module removed in order to avoid a circular module dependency,
117! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
118! added
119!
120! 1757 2016-02-22 15:49:32Z maronga
121! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
122! profiles for pressure and temperature above the LES domain.
123!
124! 1709 2015-11-04 14:47:01Z maronga
125! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
126! corrections
127!
128! 1701 2015-11-02 07:43:04Z maronga
129! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
130!
131! 1691 2015-10-26 16:17:44Z maronga
132! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
133! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
134! Added output of radiative heating rates.
135!
136! 1682 2015-10-07 23:56:08Z knoop
137! Code annotations made doxygen readable
138!
139! 1606 2015-06-29 10:43:37Z maronga
140! Added preprocessor directive __netcdf to allow for compiling without netCDF.
141! Note, however, that RRTMG cannot be used without netCDF.
142!
143! 1590 2015-05-08 13:56:27Z maronga
144! Bugfix: definition of character strings requires same length for all elements
145!
146! 1587 2015-05-04 14:19:01Z maronga
147! Added albedo class for snow
148!
149! 1585 2015-04-30 07:05:52Z maronga
150! Added support for RRTMG
151!
152! 1571 2015-03-12 16:12:49Z maronga
153! Added missing KIND attribute. Removed upper-case variable names
154!
155! 1551 2015-03-03 14:18:16Z maronga
156! Added support for data output. Various variables have been renamed. Added
157! interface for different radiation schemes (currently: clear-sky, constant, and
158! RRTM (not yet implemented).
159!
160! 1496 2014-12-02 17:25:50Z maronga
161! Initial revision
162!
163!
164! Description:
165! ------------
166!> Radiation models and interfaces
167!> @todo move variable definitions used in radiation_init only to the subroutine
168!>       as they are no longer required after initialization.
169!> @todo Output of full column vertical profiles used in RRTMG
170!> @todo Output of other rrtm arrays (such as volume mixing ratios)
171!> @todo Adapt for use with topography
172!>
173!> @note Many variables have a leading dummy dimension (0:0) in order to
174!>       match the assume-size shape expected by the RRTMG model.
175!------------------------------------------------------------------------------!
176 MODULE radiation_model_mod
177 
178    USE arrays_3d,                                                             &
179        ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw
180
181    USE cloud_parameters,                                                      &
182        ONLY:  cp, l_d_cp, rho_l
183
184    USE constants,                                                             &
185        ONLY:  pi
186
187    USE control_parameters,                                                    &
188        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
189               large_scale_forcing, lsf_surf, microphysics_morrison, phi,      & 
190               pt_surface, rho_surface, surface_pressure,                      &
191               time_since_reference_point
192
193    USE date_and_time_mod,                                                     &
194        ONLY:  calc_date_and_time, d_hours_day, d_seconds_hour, day_of_year,   &
195               time_utc
196
197    USE indices,                                                               &
198        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
199
200    USE kinds
201
202    USE microphysics_mod,                                                      &
203        ONLY:  nc_const, sigma_gc
204
205#if defined ( __netcdf )
206    USE NETCDF
207#endif
208
209#if defined ( __rrtmg )
210    USE parrrsw,                                                               &
211        ONLY:  naerec, nbndsw
212
213    USE parrrtm,                                                               &
214        ONLY:  nbndlw
215
216    USE rrtmg_lw_init,                                                         &
217        ONLY:  rrtmg_lw_ini
218
219    USE rrtmg_sw_init,                                                         &
220        ONLY:  rrtmg_sw_ini
221
222    USE rrtmg_lw_rad,                                                          &
223        ONLY:  rrtmg_lw
224
225    USE rrtmg_sw_rad,                                                          &
226        ONLY:  rrtmg_sw
227#endif
228    USE surface_mod,                                                           &
229        ONLY:  get_topography_top_index
230
231    IMPLICIT NONE
232
233    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
234
235!
236!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
237    CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
238                                   'user defined                         ', & !  0
239                                   'ocean                                ', & !  1
240                                   'mixed farming, tall grassland        ', & !  2
241                                   'tall/medium grassland                ', & !  3
242                                   'evergreen shrubland                  ', & !  4
243                                   'short grassland/meadow/shrubland     ', & !  5
244                                   'evergreen needleleaf forest          ', & !  6
245                                   'mixed deciduous evergreen forest     ', & !  7
246                                   'deciduous forest                     ', & !  8
247                                   'tropical evergreen broadleaved forest', & !  9
248                                   'medium/tall grassland/woodland       ', & ! 10
249                                   'desert, sandy                        ', & ! 11
250                                   'desert, rocky                        ', & ! 12
251                                   'tundra                               ', & ! 13
252                                   'land ice                             ', & ! 14
253                                   'sea ice                              ', & ! 15
254                                   'snow                                 ', & ! 16
255                                   'bare soil                            ', & ! 17
256                                   'asphalt/concrete mix                 ', & ! 18
257                                   'asphalt (asphalt concrete)           ', & ! 19
258                                   'concrete (Portland concrete)         ', & ! 20
259                                   'sett                                 ', & ! 21
260                                   'paving stones                        ', & ! 22
261                                   'cobblestone                          ', & ! 23
262                                   'metal                                ', & ! 24
263                                   'wood                                 ', & ! 25
264                                   'gravel                               ', & ! 26
265                                   'fine gravel                          ', & ! 27
266                                   'pebblestone                          ', & ! 28
267                                   'woodchips                            ', & ! 29
268                                   'tartan (sports)                      ', & ! 30
269                                   'artifical turf (sports)              ', & ! 31
270                                   'clay (sports)                        ', & ! 32
271                                   'building (dummy)                     '  & ! 33
272                                                         /)
273
274
275    INTEGER(iwp) :: albedo_type  = 9999999, & !< Albedo surface type
276                    dots_rad     = 0          !< starting index for timeseries output
277
278    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
279                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
280                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
281                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
282                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
283                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
284                sw_radiation = .TRUE.,                 & !< flag parameter indicing whether shortwave radiation shall be calculated
285                sun_direction = .FALSE.                 !< flag parameter indicing whether solar direction shall be calculated
286
287
288    REAL(wp), PARAMETER :: sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
289                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
290
291    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
292                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
293                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
294                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
295                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
296                decl_1,                          & !< declination coef. 1
297                decl_2,                          & !< declination coef. 2
298                decl_3,                          & !< declination coef. 3
299                dt_radiation = 0.0_wp,           & !< radiation model timestep
300                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
301                lambda = 0.0_wp,                 & !< longitude in degrees
302                lon = 0.0_wp,                    & !< longitude in radians
303                lat = 0.0_wp,                    & !< latitude in radians
304                net_radiation = 0.0_wp,          & !< net radiation at surface
305                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
306                sky_trans,                       & !< sky transmissivity
307                time_radiation = 0.0_wp            !< time since last call of radiation code
308
309
310    REAL(wp), DIMENSION(0:0) ::  zenith,         & !< cosine of solar zenith angle
311                                 sun_dir_lat,    & !< solar directional vector in latitudes
312                                 sun_dir_lon       !< solar directional vector in longitudes
313
314    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
315                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
316                emis,                        & !< surface broadband emissitivity
317                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
318                rad_net,                     & !< net radiation at the surface
319                rad_net_av                     !< average of rad_net
320
321!
322!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
323!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
324    REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 
325                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
326                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
327                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
328                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
329                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
330                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
331                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
332                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
333                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
334                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
335                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
336                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
337                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
338                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
339                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
340                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
341                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
342                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 18
343                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 19
344                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 20
345                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 21
346                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 22
347                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 23
348                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 24
349                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 25
350                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 26
351                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 27
352                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 28
353                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 29
354                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 30
355                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 31
356                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 32
357                                   0.17_wp, 0.17_wp, 0.17_wp             & ! 33
358                                 /), (/ 3, 33 /) )
359
360    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
361                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
362                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
363                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
364                        rad_lw_hr_av,                  & !< average of rad_sw_hr
365                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
366                        rad_lw_in_av,                  & !< average of rad_lw_in
367                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
368                        rad_lw_out_av,                 & !< average of rad_lw_out
369                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
370                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
371                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
372                        rad_sw_hr_av,                  & !< average of rad_sw_hr
373                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
374                        rad_sw_in_av,                  & !< average of rad_sw_in
375                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
376                        rad_sw_out_av                    !< average of rad_sw_out
377
378
379!
380!-- Variables and parameters used in RRTMG only
381#if defined ( __rrtmg )
382    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
383
384
385!
386!-- Flag parameters for RRTMGS (should not be changed)
387    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
388                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
389                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
390                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
391                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
392                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
393                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
394
395!
396!-- The following variables should be only changed with care, as this will
397!-- require further setting of some variables, which is currently not
398!-- implemented (aerosols, ice phase).
399    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
400                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
401                    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)
402
403    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
404
405    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
406
407    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
408
409    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
410                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
411                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
412                                           t_snd          !< actual temperature from sounding data (hPa)
413
414    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
415                                             aldir,          & !< longwave direct albedo solar angle of 60°
416                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
417                                             asdir,          & !< shortwave direct albedo solar angle of 60°
418                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
419                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
420                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
421                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
422                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
423                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
424                                             rrtm_cldfr,     & !< cloud fraction (0,1)
425                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
426                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
427                                             rrtm_emis,      & !< surface emissivity (0-1)   
428                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
429                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
430                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
431                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
432                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
433                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
434                                             rrtm_reice,     & !< cloud ice effective radius (microns)
435                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
436                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
437                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
438                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
439                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
440                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
441                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
442                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
443                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
444                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
445                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
446                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
447                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
448                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
449                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
450                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
451                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
452
453!
454!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
455    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
456                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
457                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
458                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
459                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
460                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
461                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
462                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
463                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
464                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
465                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
466                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
467                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
468                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
469                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
470                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
471                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
472                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
473
474#endif
475
476    INTERFACE radiation_check_data_output
477       MODULE PROCEDURE radiation_check_data_output
478    END INTERFACE radiation_check_data_output
479
480    INTERFACE radiation_check_data_output_pr
481       MODULE PROCEDURE radiation_check_data_output_pr
482    END INTERFACE radiation_check_data_output_pr
483 
484    INTERFACE radiation_check_parameters
485       MODULE PROCEDURE radiation_check_parameters
486    END INTERFACE radiation_check_parameters
487 
488    INTERFACE radiation_clearsky
489       MODULE PROCEDURE radiation_clearsky
490    END INTERFACE radiation_clearsky
491 
492    INTERFACE radiation_constant
493       MODULE PROCEDURE radiation_constant
494    END INTERFACE radiation_constant
495 
496    INTERFACE radiation_control
497       MODULE PROCEDURE radiation_control
498    END INTERFACE radiation_control
499
500    INTERFACE radiation_3d_data_averaging
501       MODULE PROCEDURE radiation_3d_data_averaging
502    END INTERFACE radiation_3d_data_averaging
503
504    INTERFACE radiation_data_output_2d
505       MODULE PROCEDURE radiation_data_output_2d
506    END INTERFACE radiation_data_output_2d
507
508    INTERFACE radiation_data_output_3d
509       MODULE PROCEDURE radiation_data_output_3d
510    END INTERFACE radiation_data_output_3d
511
512    INTERFACE radiation_data_output_mask
513       MODULE PROCEDURE radiation_data_output_mask
514    END INTERFACE radiation_data_output_mask
515
516    INTERFACE radiation_define_netcdf_grid
517       MODULE PROCEDURE radiation_define_netcdf_grid
518    END INTERFACE radiation_define_netcdf_grid
519
520    INTERFACE radiation_header
521       MODULE PROCEDURE radiation_header
522    END INTERFACE radiation_header 
523 
524    INTERFACE radiation_init
525       MODULE PROCEDURE radiation_init
526    END INTERFACE radiation_init
527
528    INTERFACE radiation_parin
529       MODULE PROCEDURE radiation_parin
530    END INTERFACE radiation_parin
531   
532    INTERFACE radiation_rrtmg
533       MODULE PROCEDURE radiation_rrtmg
534    END INTERFACE radiation_rrtmg
535
536    INTERFACE radiation_tendency
537       MODULE PROCEDURE radiation_tendency
538       MODULE PROCEDURE radiation_tendency_ij
539    END INTERFACE radiation_tendency
540
541    INTERFACE radiation_read_restart_data
542       MODULE PROCEDURE radiation_read_restart_data
543    END INTERFACE radiation_read_restart_data
544
545    INTERFACE radiation_last_actions
546       MODULE PROCEDURE radiation_last_actions
547    END INTERFACE radiation_last_actions
548
549    SAVE
550
551    PRIVATE
552
553!
554!-- Public functions / NEEDS SORTING
555    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
556           radiation_check_parameters, radiation_control,                      &
557           radiation_header, radiation_init, radiation_parin,                  &
558           radiation_3d_data_averaging, radiation_tendency,                    &
559           radiation_data_output_2d, radiation_data_output_3d,                 &
560           radiation_define_netcdf_grid, radiation_last_actions,               &
561           radiation_read_restart_data, radiation_data_output_mask
562   
563!
564!-- Public variables and constants / NEEDS SORTING
565    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation, emissivity, force_radiation_call,&
566           lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
567           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
568           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
569           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
570           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
571           skip_time_do_radiation, solar_constant, time_radiation,             &
572           unscheduled_radiation_calls, zenith, calc_zenith, sun_direction,    &
573           sun_dir_lat, sun_dir_lon
574
575
576#if defined ( __rrtmg )
577    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
578#endif
579
580 CONTAINS
581
582
583!------------------------------------------------------------------------------!
584! Description:
585! ------------
586!> This subroutine controls the calls of the radiation schemes
587!------------------------------------------------------------------------------!
588    SUBROUTINE radiation_control
589 
590 
591       IMPLICIT NONE
592
593
594       SELECT CASE ( TRIM( radiation_scheme ) )
595
596          CASE ( 'constant' )
597             CALL radiation_constant
598         
599          CASE ( 'clear-sky' )
600             CALL radiation_clearsky
601       
602          CASE ( 'rrtmg' )
603             CALL radiation_rrtmg
604
605          CASE DEFAULT
606
607       END SELECT
608
609
610    END SUBROUTINE radiation_control
611
612!------------------------------------------------------------------------------!
613! Description:
614! ------------
615!> Check data output for radiation model
616!------------------------------------------------------------------------------!
617    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
618 
619 
620       USE control_parameters,                                                 &
621           ONLY: data_output, message_string
622
623       IMPLICIT NONE
624
625       CHARACTER (LEN=*) ::  unit     !<
626       CHARACTER (LEN=*) ::  var      !<
627
628       INTEGER(iwp) :: i
629       INTEGER(iwp) :: ilen
630       INTEGER(iwp) :: k
631
632       SELECT CASE ( TRIM( var ) )
633
634         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
635             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
636                message_string = '"output of "' // TRIM( var ) // '" requi' // &
637                                 'res radiation = .TRUE. and ' //              &
638                                 'radiation_scheme = "rrtmg"'
639                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
640             ENDIF
641             unit = 'K/h'     
642
643         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
644             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
645                message_string = '"output of "' // TRIM( var ) // '" requi' // &
646                                 'res radiation = .TRUE. and ' //              &
647                                 'radiation_scheme = "rrtmg"'
648                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
649             ENDIF
650             unit = 'W/m2'   
651
652          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
653                 'rrtm_asdir*' )
654             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
655                message_string = 'illegal value for data_output: "' //         &
656                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
657                                 'cross sections are allowed for this value'
658                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
659             ENDIF
660             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
661                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
662                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
663                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
664                     TRIM( var ) == 'rrtm_asdir*'      )                       &
665                THEN
666                   message_string = 'output of "' // TRIM( var ) // '" require'&
667                                    // 's radiation = .TRUE. and radiation_sch'&
668                                    // 'eme = "rrtmg"'
669                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
670                ENDIF
671             ENDIF
672
673             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
674             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
675             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
676             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
677             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
678
679          CASE DEFAULT
680             unit = 'illegal'
681
682       END SELECT
683
684
685    END SUBROUTINE radiation_check_data_output
686
687!------------------------------------------------------------------------------!
688! Description:
689! ------------
690!> Check data output of profiles for radiation model
691!------------------------------------------------------------------------------! 
692    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
693               dopr_unit )
694 
695       USE arrays_3d,                                                          &
696           ONLY: zu
697
698       USE control_parameters,                                                 &
699           ONLY: data_output_pr, message_string
700
701       USE indices
702
703       USE profil_parameter
704
705       USE statistics
706
707       IMPLICIT NONE
708   
709       CHARACTER (LEN=*) ::  unit      !<
710       CHARACTER (LEN=*) ::  variable  !<
711       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
712 
713       INTEGER(iwp) ::  user_pr_index !<
714       INTEGER(iwp) ::  var_count     !<
715
716       SELECT CASE ( TRIM( variable ) )
717       
718         CASE ( 'rad_net' )
719             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
720             THEN
721                message_string = 'data_output_pr = ' //                        &
722                                 TRIM( data_output_pr(var_count) ) // ' is' // &
723                                 'not available for radiation = .FALSE. or ' //&
724                                 'radiation_scheme = "constant"'
725                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
726             ELSE
727                dopr_index(var_count) = 99
728                dopr_unit  = 'W/m2'
729                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
730                unit = dopr_unit
731             ENDIF
732
733          CASE ( 'rad_lw_in' )
734             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
735             THEN
736                message_string = 'data_output_pr = ' //                        &
737                                 TRIM( data_output_pr(var_count) ) // ' is' // &
738                                 'not available for radiation = .FALSE. or ' //&
739                                 'radiation_scheme = "constant"'
740                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
741             ELSE
742                dopr_index(var_count) = 100
743                dopr_unit  = 'W/m2'
744                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
745                unit = dopr_unit 
746             ENDIF
747
748          CASE ( 'rad_lw_out' )
749             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
750             THEN
751                message_string = 'data_output_pr = ' //                        &
752                                 TRIM( data_output_pr(var_count) ) // ' is' // &
753                                 'not available for radiation = .FALSE. or ' //&
754                                 'radiation_scheme = "constant"'
755                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
756             ELSE
757                dopr_index(var_count) = 101
758                dopr_unit  = 'W/m2'
759                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
760                unit = dopr_unit   
761             ENDIF
762
763          CASE ( 'rad_sw_in' )
764             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
765             THEN
766                message_string = 'data_output_pr = ' //                        &
767                                 TRIM( data_output_pr(var_count) ) // ' is' // &
768                                 'not available for radiation = .FALSE. or ' //&
769                                 'radiation_scheme = "constant"'
770                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
771             ELSE
772                dopr_index(var_count) = 102
773                dopr_unit  = 'W/m2'
774                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
775                unit = dopr_unit
776             ENDIF
777
778          CASE ( 'rad_sw_out')
779             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
780             THEN
781                message_string = 'data_output_pr = ' //                        &
782                                 TRIM( data_output_pr(var_count) ) // ' is' // &
783                                 'not available for radiation = .FALSE. or ' //&
784                                 'radiation_scheme = "constant"'
785                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
786             ELSE
787                dopr_index(var_count) = 103
788                dopr_unit  = 'W/m2'
789                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
790                unit = dopr_unit
791             ENDIF
792
793          CASE ( 'rad_lw_cs_hr' )
794             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
795             THEN
796                message_string = 'data_output_pr = ' //                        &
797                                 TRIM( data_output_pr(var_count) ) // ' is' // &
798                                 'not available for radiation = .FALSE. or ' //&
799                                 'radiation_scheme /= "rrtmg"'
800                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
801             ELSE
802                dopr_index(var_count) = 104
803                dopr_unit  = 'K/h'
804                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
805                unit = dopr_unit
806             ENDIF
807
808          CASE ( 'rad_lw_hr' )
809             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
810             THEN
811                message_string = 'data_output_pr = ' //                        &
812                                 TRIM( data_output_pr(var_count) ) // ' is' // &
813                                 'not available for radiation = .FALSE. or ' //&
814                                 'radiation_scheme /= "rrtmg"'
815                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
816             ELSE
817                dopr_index(var_count) = 105
818                dopr_unit  = 'K/h'
819                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
820                unit = dopr_unit
821             ENDIF
822
823          CASE ( 'rad_sw_cs_hr' )
824             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
825             THEN
826                message_string = 'data_output_pr = ' //                        &
827                                 TRIM( data_output_pr(var_count) ) // ' is' // &
828                                 'not available for radiation = .FALSE. or ' //&
829                                 'radiation_scheme /= "rrtmg"'
830                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
831             ELSE
832                dopr_index(var_count) = 106
833                dopr_unit  = 'K/h'
834                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
835                unit = dopr_unit
836             ENDIF
837
838          CASE ( 'rad_sw_hr' )
839             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
840             THEN
841                message_string = 'data_output_pr = ' //                        &
842                                 TRIM( data_output_pr(var_count) ) // ' is' // &
843                                 'not available for radiation = .FALSE. or ' //&
844                                 'radiation_scheme /= "rrtmg"'
845                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
846             ELSE
847                dopr_index(var_count) = 107
848                dopr_unit  = 'K/h'
849                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
850                unit = dopr_unit
851             ENDIF
852
853
854          CASE DEFAULT
855             unit = 'illegal'
856
857       END SELECT
858
859
860    END SUBROUTINE radiation_check_data_output_pr
861 
862 
863!------------------------------------------------------------------------------!
864! Description:
865! ------------
866!> Check parameters routine for radiation model
867!------------------------------------------------------------------------------!
868    SUBROUTINE radiation_check_parameters
869
870       USE control_parameters,                                                 &
871           ONLY: message_string, topography, urban_surface
872                 
873   
874       IMPLICIT NONE
875       
876
877       IF ( radiation_scheme /= 'constant'   .AND.                             &
878            radiation_scheme /= 'clear-sky'  .AND.                             &
879            radiation_scheme /= 'rrtmg' )  THEN
880          message_string = 'unknown radiation_scheme = '//                     &
881                           TRIM( radiation_scheme )
882          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
883       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
884#if ! defined ( __rrtmg )
885          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
886                           'compilation of PALM with pre-processor ' //        &
887                           'directive -D__rrtmg'
888          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
889#endif
890#if defined ( __rrtmg ) && ! defined( __netcdf )
891          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
892                           'the use of NetCDF (preprocessor directive ' //     &
893                           '-D__netcdf'
894          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
895#endif
896
897       ENDIF
898
899       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
900            radiation_scheme == 'clear-sky')  THEN
901          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
902                           'with albedo_type = 0 requires setting of albedo'// &
903                           ' /= 9999999.9'
904          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
905       ENDIF
906
907       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
908          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
909          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
910          ) ) THEN
911          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
912                           'with albedo_type = 0 requires setting of ' //      &
913                           'albedo_lw_dif /= 9999999.9' //                     &
914                           'albedo_lw_dir /= 9999999.9' //                     &
915                           'albedo_sw_dif /= 9999999.9 and' //                 &
916                           'albedo_sw_dir /= 9999999.9'
917          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
918       ENDIF
919
920!
921!--    The following paramter check is temporarily extended by the urban_surface
922!--    flag, until a better solution comes up to omit this check in case of
923!--    urban surface model is used.
924       IF ( topography /= 'flat'  .AND.  .NOT.  urban_surface )  THEN
925          message_string = 'radiation scheme cannot be used ' //               & 
926                           'in combination with  topography /= "flat"'
927          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
928       ENDIF
929 
930    END SUBROUTINE radiation_check_parameters 
931 
932 
933!------------------------------------------------------------------------------!
934! Description:
935! ------------
936!> Initialization of the radiation model
937!------------------------------------------------------------------------------!
938    SUBROUTINE radiation_init
939   
940       IMPLICIT NONE
941
942!
943!--    Allocate array for storing emissivity
944       IF ( .NOT. ALLOCATED ( emis ) )  THEN
945          ALLOCATE ( emis(nysg:nyng,nxlg:nxrg) )
946          emis = emissivity
947       ENDIF
948
949!
950!--    Allocate array for storing the surface net radiation
951       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
952          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
953          rad_net = 0.0_wp
954       ENDIF
955
956!
957!--    Allocate array for storing the surface net radiation
958       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
959          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
960          rad_lw_out_change_0 = 0.0_wp
961       ENDIF
962
963!
964!--    Fix net radiation in case of radiation_scheme = 'constant'
965       IF ( radiation_scheme == 'constant' )  THEN
966          rad_net = net_radiation
967!          radiation = .FALSE.
968!
969!--    Calculate orbital constants
970       ELSE
971          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
972          decl_2 = 2.0_wp * pi / 365.0_wp
973          decl_3 = decl_2 * 81.0_wp
974          lat    = phi * pi / 180.0_wp
975          lon    = lambda * pi / 180.0_wp
976       ENDIF
977
978       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
979            radiation_scheme == 'constant')  THEN
980
981          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
982
983          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
984             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
985          ENDIF
986          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
987             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
988          ENDIF
989
990          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
991             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
992          ENDIF
993          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
994             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
995          ENDIF
996
997          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
998             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
999          ENDIF
1000          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
1001             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
1002          ENDIF
1003
1004          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
1005             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
1006          ENDIF
1007          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
1008             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
1009          ENDIF
1010
1011          rad_sw_in  = 0.0_wp
1012          rad_sw_out = 0.0_wp
1013          rad_lw_in  = 0.0_wp
1014          rad_lw_out = 0.0_wp
1015
1016!
1017!--       Overwrite albedo if manually set in parameter file
1018          IF ( albedo_type /= 0  .AND.  albedo_type /= 9999999 .AND. albedo == 9999999.9_wp )  THEN
1019             albedo = albedo_pars(2,albedo_type)
1020          ENDIF
1021!
1022!--       Write albedo to 2d array alpha to allow surface heterogeneities   
1023          alpha = albedo
1024 
1025!
1026!--    Initialization actions for RRTMG
1027       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
1028#if defined ( __rrtmg )
1029!
1030!--       Allocate albedos
1031          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
1032          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
1033          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
1034          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
1035          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
1036          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
1037          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
1038          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
1039
1040          IF ( albedo_type /= 0 )  THEN
1041             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
1042                albedo_lw_dif = albedo_pars(0,albedo_type)
1043                albedo_lw_dir = albedo_lw_dif
1044             ENDIF
1045             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
1046                albedo_sw_dif = albedo_pars(1,albedo_type)
1047                albedo_sw_dir = albedo_sw_dif
1048             ENDIF
1049          ENDIF
1050
1051          aldif(:,:) = albedo_lw_dif
1052          aldir(:,:) = albedo_lw_dir
1053          asdif(:,:) = albedo_sw_dif
1054          asdir(:,:) = albedo_sw_dir
1055!
1056!--       Calculate initial values of current (cosine of) the zenith angle and
1057!--       whether the sun is up
1058          CALL calc_zenith     
1059!
1060!--       Calculate initial surface albedo
1061          IF ( .NOT. constant_albedo )  THEN
1062             CALL calc_albedo
1063          ELSE
1064             rrtm_aldif(0,:,:) = aldif(:,:)
1065             rrtm_aldir(0,:,:) = aldir(:,:)
1066             rrtm_asdif(0,:,:) = asdif(:,:) 
1067             rrtm_asdir(0,:,:) = asdir(:,:)   
1068          ENDIF
1069
1070!
1071!--       Allocate surface emissivity
1072          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
1073          rrtm_emis = emissivity
1074
1075!
1076!--       Allocate 3d arrays of radiative fluxes and heating rates
1077          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
1078             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1079             rad_sw_in = 0.0_wp
1080          ENDIF
1081
1082          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
1083             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1084          ENDIF
1085
1086          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
1087             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1088             rad_sw_out = 0.0_wp
1089          ENDIF
1090
1091          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
1092             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1093          ENDIF
1094
1095          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
1096             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1097             rad_sw_hr = 0.0_wp
1098          ENDIF
1099
1100          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
1101             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1102             rad_sw_hr_av = 0.0_wp
1103          ENDIF
1104
1105          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
1106             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1107             rad_sw_cs_hr = 0.0_wp
1108          ENDIF
1109
1110          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
1111             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1112             rad_sw_cs_hr_av = 0.0_wp
1113          ENDIF
1114
1115          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
1116             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1117             rad_lw_in     = 0.0_wp
1118          ENDIF
1119
1120          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
1121             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1122          ENDIF
1123
1124          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
1125             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1126            rad_lw_out    = 0.0_wp
1127          ENDIF
1128
1129          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
1130             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1131          ENDIF
1132
1133          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
1134             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1135             rad_lw_hr = 0.0_wp
1136          ENDIF
1137
1138          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
1139             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1140             rad_lw_hr_av = 0.0_wp
1141          ENDIF
1142
1143          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
1144             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1145             rad_lw_cs_hr = 0.0_wp
1146          ENDIF
1147
1148          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
1149             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1150             rad_lw_cs_hr_av = 0.0_wp
1151          ENDIF
1152
1153          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1154          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1155          rad_sw_cs_in  = 0.0_wp
1156          rad_sw_cs_out = 0.0_wp
1157
1158          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1159          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1160          rad_lw_cs_in  = 0.0_wp
1161          rad_lw_cs_out = 0.0_wp
1162
1163!
1164!--       Allocate dummy array for storing surface temperature
1165          ALLOCATE ( rrtm_tsfc(1) )
1166
1167!
1168!--       Initialize RRTMG
1169          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
1170          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
1171
1172!
1173!--       Set input files for RRTMG
1174          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
1175          IF ( .NOT. snd_exists )  THEN
1176             rrtm_input_file = "rrtmg_lw.nc"
1177          ENDIF
1178
1179!
1180!--       Read vertical layers for RRTMG from sounding data
1181!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
1182!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
1183!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
1184          CALL read_sounding_data
1185
1186!
1187!--       Read trace gas profiles from file. This routine provides
1188!--       the rrtm_ arrays (1:nzt_rad+1)
1189          CALL read_trace_gas_data
1190#endif
1191       ENDIF
1192
1193!
1194!--    Perform user actions if required
1195       CALL user_init_radiation
1196
1197!
1198!--    Calculate radiative fluxes at model start
1199       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1200
1201          SELECT CASE ( radiation_scheme )
1202             CASE ( 'rrtmg' )
1203                CALL radiation_rrtmg
1204             CASE ( 'clear-sky' )
1205                CALL radiation_clearsky
1206             CASE ( 'constant' )
1207                CALL radiation_constant
1208             CASE DEFAULT
1209          END SELECT
1210
1211       ENDIF
1212
1213       RETURN
1214
1215    END SUBROUTINE radiation_init
1216
1217
1218!------------------------------------------------------------------------------!
1219! Description:
1220! ------------
1221!> A simple clear sky radiation model
1222!------------------------------------------------------------------------------!
1223    SUBROUTINE radiation_clearsky
1224
1225
1226       IMPLICIT NONE
1227
1228       INTEGER(iwp) :: i, j, k   !< loop indices
1229       REAL(wp)     :: exn,   &  !< Exner functions at surface
1230                       exn1,  &  !< Exner functions at first grid level
1231                       pt1       !< potential temperature at first grid level
1232
1233!
1234!--    Calculate current zenith angle
1235       CALL calc_zenith
1236
1237!
1238!--    Calculate sky transmissivity
1239       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1240
1241!
1242!--    Calculate value of the Exner function
1243       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1244!
1245!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
1246!--    point
1247       DO i = nxlg, nxrg
1248          DO j = nysg, nyng
1249!
1250!--          Obtain vertical index of topography top
1251             k = get_topography_top_index( j, i, 's' )
1252
1253             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1254
1255             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
1256             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1257             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
1258
1259             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
1260                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1261                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1262             ELSE
1263                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1264             ENDIF
1265
1266             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
1267                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1268
1269
1270             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emis(j,i)         &
1271                                        * (pt(k,j,i) * exn) ** 3
1272
1273          ENDDO
1274       ENDDO
1275
1276    END SUBROUTINE radiation_clearsky
1277
1278
1279!------------------------------------------------------------------------------!
1280! Description:
1281! ------------
1282!> This scheme keeps the prescribed net radiation constant during the run
1283!------------------------------------------------------------------------------!
1284    SUBROUTINE radiation_constant
1285
1286
1287       IMPLICIT NONE
1288
1289       INTEGER(iwp) :: i, j, k   !< loop indices
1290       REAL(wp)     :: exn,   &  !< Exner functions at surface
1291                       exn1,  &  !< Exner functions at first grid level
1292                       pt1       !< potential temperature at first grid level
1293
1294!
1295!--    Calculate value of the Exner function
1296       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1297!
1298!--    Prescribe net radiation and estimate the remaining radiative fluxes
1299       DO i = nxlg, nxrg
1300          DO j = nysg, nyng
1301!
1302!--          Obtain vertical index of topography top. So far it is identical to
1303!--          nzb.
1304             k = get_topography_top_index( j, i, 's' )
1305
1306             rad_net(j,i)      = net_radiation
1307
1308             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1309
1310             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
1311                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1312                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1313             ELSE
1314                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1315             ENDIF
1316
1317             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
1318
1319             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
1320                                  + rad_lw_out(0,j,i) )                        &
1321                                  / ( 1.0_wp - alpha(j,i) )
1322
1323             rad_sw_out(0,j,i) =  alpha(j,i) * rad_sw_in(0,j,i)
1324
1325          ENDDO
1326       ENDDO
1327
1328    END SUBROUTINE radiation_constant
1329
1330!------------------------------------------------------------------------------!
1331! Description:
1332! ------------
1333!> Header output for radiation model
1334!------------------------------------------------------------------------------!
1335    SUBROUTINE radiation_header ( io )
1336
1337
1338       IMPLICIT NONE
1339 
1340       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1341   
1342
1343       
1344!
1345!--    Write radiation model header
1346       WRITE( io, 3 )
1347
1348       IF ( radiation_scheme == "constant" )  THEN
1349          WRITE( io, 4 ) net_radiation
1350       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1351          WRITE( io, 5 )
1352       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1353          WRITE( io, 6 )
1354          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1355          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1356       ENDIF
1357
1358       IF ( albedo_type == 0 )  THEN
1359          WRITE( io, 7 ) albedo
1360       ELSE
1361          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1362       ENDIF
1363       IF ( constant_albedo )  THEN
1364          WRITE( io, 9 )
1365       ENDIF
1366       
1367       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1368          WRITE ( io, 1 )  lambda
1369       ENDIF
1370
1371       WRITE( io, 12 ) dt_radiation
1372 
1373
1374 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1375 3 FORMAT (//' Radiation model information:'/                                  &
1376              ' ----------------------------'/)
1377 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
1378           // 'W/m**2')
1379 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
1380                   ' default)')
1381 6 FORMAT ('    --> RRTMG scheme is used')
1382 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1383 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1384 9 FORMAT (/'    --> Albedo is fixed during the run')
138510 FORMAT (/'    --> Longwave radiation is disabled')
138611 FORMAT (/'    --> Shortwave radiation is disabled.')
138712 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1388
1389
1390    END SUBROUTINE radiation_header
1391   
1392
1393!------------------------------------------------------------------------------!
1394! Description:
1395! ------------
1396!> Parin for &radiation_par for radiation model
1397!------------------------------------------------------------------------------!
1398    SUBROUTINE radiation_parin
1399
1400
1401       IMPLICIT NONE
1402
1403       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1404       
1405       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
1406                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1407                                  constant_albedo, dt_radiation,               &
1408                                  lambda, lw_radiation, net_radiation,         &
1409                                  radiation_scheme, skip_time_do_radiation,    &
1410                                  sw_radiation, unscheduled_radiation_calls
1411       
1412       line = ' '
1413       
1414!
1415!--    Try to find radiation model package
1416       REWIND ( 11 )
1417       line = ' '
1418       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1419          READ ( 11, '(A)', END=10 )  line
1420       ENDDO
1421       BACKSPACE ( 11 )
1422
1423!
1424!--    Read user-defined namelist
1425       READ ( 11, radiation_par )
1426
1427!
1428!--    Set flag that indicates that the radiation model is switched on
1429       radiation = .TRUE.
1430
1431 10    CONTINUE
1432       
1433
1434    END SUBROUTINE radiation_parin
1435
1436
1437!------------------------------------------------------------------------------!
1438! Description:
1439! ------------
1440!> Implementation of the RRTMG radiation_scheme
1441!------------------------------------------------------------------------------!
1442    SUBROUTINE radiation_rrtmg
1443
1444       USE indices,                                                            &
1445           ONLY:  nbgp
1446
1447       USE particle_attributes,                                                &
1448           ONLY:  grid_particles, number_of_particles, particles,              &
1449                  particle_advection_start, prt_count
1450
1451       IMPLICIT NONE
1452
1453#if defined ( __rrtmg )
1454
1455       INTEGER(iwp) :: i, j, k, n !< loop indices
1456
1457       REAL(wp)     ::  nc_rad, &    !< number concentration of cloud droplets
1458                        s_r2,   &    !< weighted sum over all droplets with r^2
1459                        s_r3         !< weighted sum over all droplets with r^3
1460
1461!
1462!--    Calculate current (cosine of) zenith angle and whether the sun is up
1463       CALL calc_zenith     
1464!
1465!--    Calculate surface albedo
1466       IF ( .NOT. constant_albedo )  THEN
1467          CALL calc_albedo
1468       ENDIF
1469
1470!
1471!--    Prepare input data for RRTMG
1472
1473!
1474!--    In case of large scale forcing with surface data, calculate new pressure
1475!--    profile. nzt_rad might be modified by these calls and all required arrays
1476!--    will then be re-allocated
1477       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1478          CALL read_sounding_data
1479          CALL read_trace_gas_data
1480       ENDIF
1481!
1482!--    Loop over all grid points
1483       DO i = nxl, nxr
1484          DO j = nys, nyn
1485
1486!
1487!--          Prepare profiles of temperature and H2O volume mixing ratio
1488             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
1489                                                  / 1000.0_wp )**0.286_wp
1490
1491
1492             IF ( cloud_physics )  THEN
1493                DO k = nzb+1, nzt+1
1494                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1495                                    )**0.286_wp + l_d_cp * ql(k,j,i)
1496                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1497                ENDDO
1498             ELSEIF ( cloud_droplets ) THEN
1499                DO k = nzb+1, nzt+1
1500                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1501                                    )**0.286_wp + l_d_cp * ql(k,j,i)
1502                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 
1503                ENDDO     
1504             ELSE
1505                DO k = nzb+1, nzt+1
1506                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1507                                    )**0.286_wp
1508                   rrtm_h2ovmr(0,k) = 0.0_wp
1509                ENDDO
1510             ENDIF
1511
1512!
1513!--          Avoid temperature/humidity jumps at the top of the LES domain by
1514!--          linear interpolation from nzt+2 to nzt+7
1515             DO k = nzt+2, nzt+7
1516                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
1517                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
1518                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
1519                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1520
1521                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
1522                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1523                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
1524                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1525
1526             ENDDO
1527
1528!--          Linear interpolate to zw grid
1529             DO k = nzb+2, nzt+8
1530                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
1531                                   rrtm_tlay(0,k-1))                           &
1532                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
1533                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1534             ENDDO
1535
1536
1537!
1538!--          Calculate liquid water path and cloud fraction for each column.
1539!--          Note that LWP is required in g/m² instead of kg/kg m.
1540             rrtm_cldfr  = 0.0_wp
1541             rrtm_reliq  = 0.0_wp
1542             rrtm_cliqwp = 0.0_wp
1543             rrtm_icld   = 0
1544
1545             IF ( cloud_physics  .OR.  cloud_droplets )  THEN
1546                DO k = nzb+1, nzt+1
1547                   rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                 &
1548                                       (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
1549                                       * 100.0_wp / g 
1550
1551                   IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
1552                      rrtm_cldfr(0,k) = 1.0_wp
1553                      IF ( rrtm_icld == 0 )  rrtm_icld = 1
1554
1555!
1556!--                   Calculate cloud droplet effective radius
1557                      IF ( cloud_physics )  THEN
1558                         nc_rad = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
1559                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
1560                                           * rho_surface                       &
1561                                           / ( 4.0_wp * pi * nc_rad * rho_l )&
1562                                           )**0.33333333333333_wp              &
1563                                           * EXP( LOG( sigma_gc )**2 )
1564
1565                      ELSEIF ( cloud_droplets )  THEN
1566                         number_of_particles = prt_count(k,j,i)
1567
1568                         IF (number_of_particles <= 0)  CYCLE
1569                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1570                         s_r2 = 0.0_wp
1571                         s_r3 = 0.0_wp
1572
1573                         DO  n = 1, number_of_particles
1574                            IF ( particles(n)%particle_mask )  THEN
1575                               s_r2 = s_r2 + particles(n)%radius**2 * &
1576                                      particles(n)%weight_factor
1577                               s_r3 = s_r3 + particles(n)%radius**3 * &
1578                                      particles(n)%weight_factor
1579                            ENDIF
1580                         ENDDO
1581
1582                         IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
1583
1584                      ENDIF
1585
1586!
1587!--                   Limit effective radius
1588                      IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
1589                         rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1590                         rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1591                     ENDIF
1592                   ENDIF
1593                ENDDO
1594             ENDIF
1595
1596!
1597!--          Set surface temperature
1598             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1599
1600!
1601!--          Set surface emissivity
1602             rrtm_emis = emis(j,i)
1603
1604             IF ( lw_radiation )  THEN
1605               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
1606               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
1607               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
1608               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
1609               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
1610               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
1611               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
1612               rrtm_reliq      , rrtm_lw_tauaer,                               &
1613               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
1614               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
1615               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
1616
1617!
1618!--             Save fluxes
1619                DO k = nzb, nzt+1
1620                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
1621                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1622                ENDDO
1623
1624!
1625!--             Save heating rates (convert from K/d to K/h)
1626                DO k = nzb+1, nzt+1
1627                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
1628                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
1629                ENDDO
1630
1631!
1632!--             Save change in LW heating rate
1633                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1634
1635             ENDIF
1636
1637             IF ( sw_radiation .AND. sun_up )  THEN
1638                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
1639               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
1640               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
1641               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1642               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1643               0.0_wp          , day_of_year  , solar_constant,   rrtm_inflgsw,&
1644               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
1645               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
1646               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
1647               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
1648               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
1649               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
1650 
1651!
1652!--             Save fluxes
1653                DO k = nzb, nzt+1
1654                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
1655                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1656                ENDDO
1657
1658!
1659!--             Save heating rates (convert from K/d to K/s)
1660                DO k = nzb+1, nzt+1
1661                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
1662                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
1663                ENDDO
1664
1665             ENDIF
1666
1667!
1668!--          Calculate surface net radiation
1669             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
1670                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1671
1672          ENDDO
1673       ENDDO
1674
1675       CALL exchange_horiz( rad_lw_in,  nbgp )
1676       CALL exchange_horiz( rad_lw_out, nbgp )
1677       CALL exchange_horiz( rad_lw_hr,    nbgp )
1678       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1679
1680       CALL exchange_horiz( rad_sw_in,  nbgp )
1681       CALL exchange_horiz( rad_sw_out, nbgp ) 
1682       CALL exchange_horiz( rad_sw_hr,    nbgp )
1683       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1684
1685       CALL exchange_horiz_2d( rad_net )
1686       CALL exchange_horiz_2d( rad_lw_out_change_0 )
1687#endif
1688
1689    END SUBROUTINE radiation_rrtmg
1690
1691
1692!------------------------------------------------------------------------------!
1693! Description:
1694! ------------
1695!> Calculate the cosine of the zenith angle (variable is called zenith)
1696!------------------------------------------------------------------------------!
1697    SUBROUTINE calc_zenith
1698
1699       IMPLICIT NONE
1700
1701       REAL(wp) ::  declination,  & !< solar declination angle
1702                    hour_angle      !< solar hour angle
1703!
1704!--    Calculate current day and time based on the initial values and simulation
1705!--    time
1706       CALL calc_date_and_time
1707
1708!
1709!--    Calculate solar declination and hour angle   
1710       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day_of_year, KIND=wp) - decl_3) )
1711       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1712
1713!
1714!--    Calculate cosine of solar zenith angle
1715       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
1716                                            * COS(hour_angle)
1717       zenith(0) = MAX(0.0_wp,zenith(0))
1718
1719!
1720!--    Calculate solar directional vector
1721       IF ( sun_direction )  THEN
1722
1723!
1724!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
1725          sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
1726
1727!
1728!--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
1729          sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
1730                              * COS(declination) * SIN(lat)
1731       ENDIF
1732
1733!
1734!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
1735       IF ( zenith(0) > 0.0_wp )  THEN
1736          sun_up = .TRUE.
1737       ELSE
1738          sun_up = .FALSE.
1739       END IF
1740
1741    END SUBROUTINE calc_zenith
1742
1743#if defined ( __rrtmg ) && defined ( __netcdf )
1744!------------------------------------------------------------------------------!
1745! Description:
1746! ------------
1747!> Calculates surface albedo components based on Briegleb (1992) and
1748!> Briegleb et al. (1986)
1749!------------------------------------------------------------------------------!
1750    SUBROUTINE calc_albedo
1751
1752        IMPLICIT NONE
1753
1754        IF ( sun_up )  THEN
1755!
1756!--        Ocean
1757           IF ( albedo_type == 1 )  THEN
1758              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
1759                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
1760                                            * ( zenith(0) - 0.5_wp )           &
1761                                            * ( zenith(0) - 1.0_wp )
1762              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1763!
1764!--        Snow
1765           ELSEIF ( albedo_type == 16 )  THEN
1766              IF ( zenith(0) < 0.5_wp )  THEN
1767                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
1768                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1769                                     * zenith(0))) - 1.0_wp
1770                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
1771                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1772                                     * zenith(0))) - 1.0_wp
1773
1774                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1775                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1776              ELSE
1777                 rrtm_aldir(0,:,:) = aldif
1778                 rrtm_asdir(0,:,:) = asdif
1779              ENDIF
1780!
1781!--        Sea ice
1782           ELSEIF ( albedo_type == 15 )  THEN
1783                 rrtm_aldir(0,:,:) = aldif
1784                 rrtm_asdir(0,:,:) = asdif
1785
1786!
1787!--        Bare soil
1788           ELSEIF ( albedo_type == 17 )  THEN
1789                 rrtm_aldir(0,:,:) = aldif
1790                 rrtm_asdir(0,:,:) = asdif
1791
1792!
1793!--        For impermeable surfaces, use values from the lookup table
1794           ELSEIF ( albedo_type > 17 )  THEN
1795                 rrtm_aldir(0,:,:) = aldif
1796                 rrtm_asdir(0,:,:) = asdif
1797!
1798!--        Land surfaces
1799           ELSE
1800               SELECT CASE ( albedo_type )
1801
1802!
1803!--              Surface types with strong zenith dependence
1804                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
1805                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
1806                                        (1.0_wp + 0.8_wp * zenith(0))
1807                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
1808                                        (1.0_wp + 0.8_wp * zenith(0))
1809!
1810!--              Surface types with weak zenith dependence
1811                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
1812                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
1813                                        (1.0_wp + 0.2_wp * zenith(0))
1814                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
1815                                        (1.0_wp + 0.2_wp * zenith(0))
1816
1817                 CASE DEFAULT
1818
1819              END SELECT
1820           ENDIF
1821!
1822!--        Diffusive albedo is taken from Table 2
1823           rrtm_aldif(0,:,:) = aldif
1824           rrtm_asdif(0,:,:) = asdif
1825
1826        ELSE
1827
1828           rrtm_aldir(0,:,:) = 0.0_wp
1829           rrtm_asdir(0,:,:) = 0.0_wp
1830           rrtm_aldif(0,:,:) = 0.0_wp
1831           rrtm_asdif(0,:,:) = 0.0_wp
1832        ENDIF
1833    END SUBROUTINE calc_albedo
1834
1835!------------------------------------------------------------------------------!
1836! Description:
1837! ------------
1838!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1839!------------------------------------------------------------------------------!
1840    SUBROUTINE read_sounding_data
1841
1842       IMPLICIT NONE
1843
1844       INTEGER(iwp) :: id,           & !< NetCDF id of input file
1845                       id_dim_zrad,  & !< pressure level id in the NetCDF file
1846                       id_var,       & !< NetCDF variable id
1847                       k,            & !< loop index
1848                       nz_snd,       & !< number of vertical levels in the sounding data
1849                       nz_snd_start, & !< start vertical index for sounding data to be used
1850                       nz_snd_end      !< end vertical index for souding data to be used
1851
1852       REAL(wp) :: t_surface           !< actual surface temperature
1853
1854       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1855                                               t_snd_tmp      !< temporary temperature profile (sounding)
1856
1857!
1858!--    In case of updates, deallocate arrays first (sufficient to check one
1859!--    array as the others are automatically allocated). This is required
1860!--    because nzt_rad might change during the update
1861       IF ( ALLOCATED ( hyp_snd ) )  THEN
1862          DEALLOCATE( hyp_snd )
1863          DEALLOCATE( t_snd )
1864          DEALLOCATE( q_snd  )
1865          DEALLOCATE ( rrtm_play )
1866          DEALLOCATE ( rrtm_plev )
1867          DEALLOCATE ( rrtm_tlay )
1868          DEALLOCATE ( rrtm_tlev )
1869
1870          DEALLOCATE ( rrtm_h2ovmr )
1871          DEALLOCATE ( rrtm_cicewp )
1872          DEALLOCATE ( rrtm_cldfr )
1873          DEALLOCATE ( rrtm_cliqwp )
1874          DEALLOCATE ( rrtm_reice )
1875          DEALLOCATE ( rrtm_reliq )
1876          DEALLOCATE ( rrtm_lw_taucld )
1877          DEALLOCATE ( rrtm_lw_tauaer )
1878
1879          DEALLOCATE ( rrtm_lwdflx  )
1880          DEALLOCATE ( rrtm_lwdflxc )
1881          DEALLOCATE ( rrtm_lwuflx  )
1882          DEALLOCATE ( rrtm_lwuflxc )
1883          DEALLOCATE ( rrtm_lwuflx_dt )
1884          DEALLOCATE ( rrtm_lwuflxc_dt )
1885          DEALLOCATE ( rrtm_lwhr  )
1886          DEALLOCATE ( rrtm_lwhrc )
1887
1888          DEALLOCATE ( rrtm_sw_taucld )
1889          DEALLOCATE ( rrtm_sw_ssacld )
1890          DEALLOCATE ( rrtm_sw_asmcld )
1891          DEALLOCATE ( rrtm_sw_fsfcld )
1892          DEALLOCATE ( rrtm_sw_tauaer )
1893          DEALLOCATE ( rrtm_sw_ssaaer )
1894          DEALLOCATE ( rrtm_sw_asmaer ) 
1895          DEALLOCATE ( rrtm_sw_ecaer )   
1896 
1897          DEALLOCATE ( rrtm_swdflx  )
1898          DEALLOCATE ( rrtm_swdflxc )
1899          DEALLOCATE ( rrtm_swuflx  )
1900          DEALLOCATE ( rrtm_swuflxc )
1901          DEALLOCATE ( rrtm_swhr  )
1902          DEALLOCATE ( rrtm_swhrc )
1903
1904       ENDIF
1905
1906!
1907!--    Open file for reading
1908       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1909       CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1910
1911!
1912!--    Inquire dimension of z axis and save in nz_snd
1913       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1914       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1915       CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1916
1917!
1918! !--    Allocate temporary array for storing pressure data
1919       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1920       hyp_snd_tmp = 0.0_wp
1921
1922
1923!--    Read pressure from file
1924       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1925       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
1926                               count = (/nz_snd/) )
1927       CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1928
1929!
1930!--    Allocate temporary array for storing temperature data
1931       ALLOCATE( t_snd_tmp(1:nz_snd) )
1932       t_snd_tmp = 0.0_wp
1933
1934!
1935!--    Read temperature from file
1936       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1937       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
1938                               count = (/nz_snd/) )
1939       CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1940
1941!
1942!--    Calculate start of sounding data
1943       nz_snd_start = nz_snd + 1
1944       nz_snd_end   = nz_snd + 1
1945
1946!
1947!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
1948!--    in Pa, hyp_snd in hPa).
1949       DO  k = 1, nz_snd
1950          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
1951             nz_snd_start = k
1952             EXIT
1953          END IF
1954       END DO
1955
1956       IF ( nz_snd_start <= nz_snd )  THEN
1957          nz_snd_end = nz_snd
1958       END IF
1959
1960
1961!
1962!--    Calculate of total grid points for RRTMG calculations
1963       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1964
1965!
1966!--    Save data above LES domain in hyp_snd, t_snd and q_snd
1967!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1968       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1969       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
1970       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
1971       hyp_snd = 0.0_wp
1972       t_snd = 0.0_wp
1973       q_snd = 0.0_wp
1974
1975       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1976       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1977
1978       nc_stat = NF90_CLOSE( id )
1979
1980!
1981!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
1982!--    top of the LES domain. This routine does not consider horizontal or
1983!--    vertical variability of pressure and temperature
1984       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
1985       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
1986
1987       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1988       DO k = nzb+1, nzt+1
1989          rrtm_play(0,k) = hyp(k) * 0.01_wp
1990          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1991                         t_surface )**(1.0_wp/0.286_wp)
1992       ENDDO
1993
1994       DO k = nzt+2, nzt_rad
1995          rrtm_play(0,k) = hyp_snd(k)
1996          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1997       ENDDO
1998       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
1999                                   1.5 * hyp_snd(nzt_rad)                      &
2000                                 - 0.5 * hyp_snd(nzt_rad-1) )
2001       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
2002                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
2003
2004       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
2005
2006!
2007!--    Calculate temperature/humidity levels at top of the LES domain.
2008!--    Currently, the temperature is taken from sounding data (might lead to a
2009!--    temperature jump at interface. To do: Humidity is currently not
2010!--    calculated above the LES domain.
2011       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
2012       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
2013       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
2014
2015       DO k = nzt+8, nzt_rad
2016          rrtm_tlay(0,k)   = t_snd(k)
2017          rrtm_h2ovmr(0,k) = q_snd(k)
2018       ENDDO
2019       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
2020                                - rrtm_tlay(0,nzt_rad-1)
2021       DO k = nzt+9, nzt_rad+1
2022          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
2023                             - rrtm_tlay(0,k-1))                               &
2024                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
2025                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
2026       ENDDO
2027       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
2028
2029       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
2030                                  - rrtm_tlev(0,nzt_rad)
2031!
2032!--    Allocate remaining RRTMG arrays
2033       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
2034       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
2035       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
2036       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
2037       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
2038       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
2039       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
2040       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
2041       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
2042       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
2043       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
2044       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
2045       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
2046       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
2047       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
2048
2049!
2050!--    The ice phase is currently not considered in PALM
2051       rrtm_cicewp = 0.0_wp
2052       rrtm_reice  = 0.0_wp
2053
2054!
2055!--    Set other parameters (move to NAMELIST parameters in the future)
2056       rrtm_lw_tauaer = 0.0_wp
2057       rrtm_lw_taucld = 0.0_wp
2058       rrtm_sw_taucld = 0.0_wp
2059       rrtm_sw_ssacld = 0.0_wp
2060       rrtm_sw_asmcld = 0.0_wp
2061       rrtm_sw_fsfcld = 0.0_wp
2062       rrtm_sw_tauaer = 0.0_wp
2063       rrtm_sw_ssaaer = 0.0_wp
2064       rrtm_sw_asmaer = 0.0_wp
2065       rrtm_sw_ecaer  = 0.0_wp
2066
2067
2068       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
2069       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
2070       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
2071       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
2072       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
2073       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
2074
2075       rrtm_swdflx  = 0.0_wp
2076       rrtm_swuflx  = 0.0_wp
2077       rrtm_swhr    = 0.0_wp 
2078       rrtm_swuflxc = 0.0_wp
2079       rrtm_swdflxc = 0.0_wp
2080       rrtm_swhrc   = 0.0_wp
2081
2082       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
2083       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
2084       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
2085       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
2086       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
2087       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
2088
2089       rrtm_lwdflx  = 0.0_wp
2090       rrtm_lwuflx  = 0.0_wp
2091       rrtm_lwhr    = 0.0_wp 
2092       rrtm_lwuflxc = 0.0_wp
2093       rrtm_lwdflxc = 0.0_wp
2094       rrtm_lwhrc   = 0.0_wp
2095
2096       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
2097       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
2098
2099       rrtm_lwuflx_dt = 0.0_wp
2100       rrtm_lwuflxc_dt = 0.0_wp
2101
2102    END SUBROUTINE read_sounding_data
2103
2104
2105!------------------------------------------------------------------------------!
2106! Description:
2107! ------------
2108!> Read trace gas data from file
2109!------------------------------------------------------------------------------!
2110    SUBROUTINE read_trace_gas_data
2111
2112       USE rrsw_ncpar
2113
2114       IMPLICIT NONE
2115
2116       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
2117
2118       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
2119           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
2120                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
2121
2122       INTEGER(iwp) :: id,     & !< NetCDF id
2123                       k,      & !< loop index
2124                       m,      & !< loop index
2125                       n,      & !< loop index
2126                       nabs,   & !< number of absorbers
2127                       np,     & !< number of pressure levels
2128                       id_abs, & !< NetCDF id of the respective absorber
2129                       id_dim, & !< NetCDF id of asborber's dimension
2130                       id_var    !< NetCDf id ot the absorber
2131
2132       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
2133
2134
2135       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,          & !< pressure levels for the absorbers
2136                                                 rrtm_play_tmp,  & !< temporary array for pressure zu-levels
2137                                                 rrtm_plev_tmp,  & !< temporary array for pressure zw-levels
2138                                                 trace_path_tmp    !< temporary array for storing trace gas path data
2139
2140       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
2141                                                 trace_mls_path, & !< array for storing trace gas path data
2142                                                 trace_mls_tmp     !< temporary array for storing trace gas data
2143
2144
2145!
2146!--    In case of updates, deallocate arrays first (sufficient to check one
2147!--    array as the others are automatically allocated)
2148       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
2149          DEALLOCATE ( rrtm_o3vmr  )
2150          DEALLOCATE ( rrtm_co2vmr )
2151          DEALLOCATE ( rrtm_ch4vmr )
2152          DEALLOCATE ( rrtm_n2ovmr )
2153          DEALLOCATE ( rrtm_o2vmr  )
2154          DEALLOCATE ( rrtm_cfc11vmr )
2155          DEALLOCATE ( rrtm_cfc12vmr )
2156          DEALLOCATE ( rrtm_cfc22vmr )
2157          DEALLOCATE ( rrtm_ccl4vmr  )
2158       ENDIF
2159
2160!
2161!--    Allocate trace gas profiles
2162       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
2163       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
2164       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
2165       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
2166       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
2167       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
2168       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
2169       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
2170       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
2171
2172!
2173!--    Open file for reading
2174       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
2175       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
2176!
2177!--    Inquire dimension ids and dimensions
2178       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
2179       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2180       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
2181       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2182
2183       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
2184       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2185       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
2186       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2187   
2188
2189!
2190!--    Allocate pressure, and trace gas arrays     
2191       ALLOCATE( p_mls(1:np) )
2192       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
2193       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
2194
2195
2196       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
2197       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2198       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
2199       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2200
2201       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
2202       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2203       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
2204       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2205
2206
2207!
2208!--    Write absorber amounts (mls) to trace_mls
2209       DO n = 1, num_trace_gases
2210          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
2211
2212          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
2213
2214!
2215!--       Replace missing values by zero
2216          WHERE ( trace_mls(n,:) > 2.0_wp ) 
2217             trace_mls(n,:) = 0.0_wp
2218          END WHERE
2219       END DO
2220
2221       DEALLOCATE ( trace_mls_tmp )
2222
2223       nc_stat = NF90_CLOSE( id )
2224       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
2225
2226!
2227!--    Add extra pressure level for calculations of the trace gas paths
2228       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
2229       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
2230
2231       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
2232       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
2233       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
2234       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
2235                                         * rrtm_plev(0,nzt_rad+1) )
2236 
2237!
2238!--    Calculate trace gas path (zero at surface) with interpolation to the
2239!--    sounding levels
2240       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
2241
2242       trace_mls_path(nzb+1,:) = 0.0_wp
2243       
2244       DO k = nzb+2, nzt_rad+2
2245          DO m = 1, num_trace_gases
2246             trace_mls_path(k,m) = trace_mls_path(k-1,m)
2247
2248!
2249!--          When the pressure level is higher than the trace gas pressure
2250!--          level, assume that
2251             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
2252               
2253                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
2254                                      * ( rrtm_plev_tmp(k-1)                   &
2255                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
2256                                        ) / g
2257             ENDIF
2258
2259!
2260!--          Integrate for each sounding level from the contributing p_mls
2261!--          levels
2262             DO n = 2, np
2263!
2264!--             Limit p_mls so that it is within the model level
2265                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
2266                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
2267                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
2268                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
2269
2270                IF ( p_mls_l > p_mls_u )  THEN
2271
2272!
2273!--                Calculate weights for interpolation
2274                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
2275                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
2276                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
2277
2278!
2279!--                Add level to trace gas path
2280                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
2281                                         +  ( p_wgt_u * trace_mls(m,n)         &
2282                                            + p_wgt_l * trace_mls(m,n-1) )     &
2283                                         * (p_mls_l - p_mls_u) / g
2284                ENDIF
2285             ENDDO
2286
2287             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
2288                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
2289                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
2290                                          - rrtm_plev_tmp(k)                   &
2291                                        ) / g 
2292             ENDIF 
2293          ENDDO
2294       ENDDO
2295
2296
2297!
2298!--    Prepare trace gas path profiles
2299       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
2300
2301       DO m = 1, num_trace_gases
2302
2303          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
2304                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
2305                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
2306                                       - rrtm_plev_tmp(2:nzt_rad+2) )
2307
2308!
2309!--       Save trace gas paths to the respective arrays
2310          SELECT CASE ( TRIM( trace_names(m) ) )
2311
2312             CASE ( 'O3' )
2313
2314                rrtm_o3vmr(0,:) = trace_path_tmp(:)
2315
2316             CASE ( 'CO2' )
2317
2318                rrtm_co2vmr(0,:) = trace_path_tmp(:)
2319
2320             CASE ( 'CH4' )
2321
2322                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2323
2324             CASE ( 'N2O' )
2325
2326                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2327
2328             CASE ( 'O2' )
2329
2330                rrtm_o2vmr(0,:) = trace_path_tmp(:)
2331
2332             CASE ( 'CFC11' )
2333
2334                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2335
2336             CASE ( 'CFC12' )
2337
2338                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2339
2340             CASE ( 'CFC22' )
2341
2342                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2343
2344             CASE ( 'CCL4' )
2345
2346                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2347
2348             CASE DEFAULT
2349
2350          END SELECT
2351
2352       ENDDO
2353
2354       DEALLOCATE ( trace_path_tmp )
2355       DEALLOCATE ( trace_mls_path )
2356       DEALLOCATE ( rrtm_play_tmp )
2357       DEALLOCATE ( rrtm_plev_tmp )
2358       DEALLOCATE ( trace_mls )
2359       DEALLOCATE ( p_mls )
2360
2361    END SUBROUTINE read_trace_gas_data
2362
2363
2364    SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2365
2366       USE control_parameters,                                                 &
2367           ONLY:  message_string
2368
2369       USE NETCDF
2370
2371       USE pegrid
2372
2373       IMPLICIT NONE
2374
2375       CHARACTER(LEN=6) ::  message_identifier
2376       CHARACTER(LEN=*) ::  routine_name
2377
2378       INTEGER(iwp) ::  errno
2379
2380       IF ( nc_stat /= NF90_NOERR )  THEN
2381
2382          WRITE( message_identifier, '(''NC'',I4.4)' )  errno
2383          message_string = TRIM( NF90_STRERROR( nc_stat ) )
2384
2385          CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2386
2387       ENDIF
2388
2389    END SUBROUTINE netcdf_handle_error_rad
2390#endif
2391
2392
2393!------------------------------------------------------------------------------!
2394! Description:
2395! ------------
2396!> Calculate temperature tendency due to radiative cooling/heating.
2397!> Cache-optimized version.
2398!------------------------------------------------------------------------------!
2399 SUBROUTINE radiation_tendency_ij ( i, j, tend )
2400
2401    USE cloud_parameters,                                                      &
2402        ONLY:  pt_d_t
2403
2404    IMPLICIT NONE
2405
2406    INTEGER(iwp) :: i, j, k !< loop indices
2407
2408    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2409
2410    IF ( radiation_scheme == 'rrtmg' )  THEN
2411#if defined  ( __rrtmg )
2412!
2413!--    Calculate tendency based on heating rate
2414       DO k = nzb+1, nzt+1
2415          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
2416                                         * pt_d_t(k) * d_seconds_hour
2417       ENDDO
2418#endif
2419    ENDIF
2420
2421    END SUBROUTINE radiation_tendency_ij
2422
2423
2424!------------------------------------------------------------------------------!
2425! Description:
2426! ------------
2427!> Calculate temperature tendency due to radiative cooling/heating.
2428!> Vector-optimized version
2429!------------------------------------------------------------------------------!
2430 SUBROUTINE radiation_tendency ( tend )
2431
2432    USE cloud_parameters,                                                      &
2433        ONLY:  pt_d_t
2434
2435    USE indices,                                                               &
2436        ONLY:  nxl, nxr, nyn, nys
2437
2438    IMPLICIT NONE
2439
2440    INTEGER(iwp) :: i, j, k !< loop indices
2441
2442    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2443
2444    IF ( radiation_scheme == 'rrtmg' )  THEN
2445#if defined  ( __rrtmg )
2446!
2447!--    Calculate tendency based on heating rate
2448       DO  i = nxl, nxr
2449          DO  j = nys, nyn
2450             DO k = nzb+1, nzt+1
2451                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
2452                                          +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
2453                                          * d_seconds_hour
2454             ENDDO
2455          ENDDO
2456       ENDDO
2457#endif
2458    ENDIF
2459
2460
2461 END SUBROUTINE radiation_tendency
2462
2463!------------------------------------------------------------------------------!
2464!
2465! Description:
2466! ------------
2467!> Subroutine for averaging 3D data
2468!------------------------------------------------------------------------------!
2469SUBROUTINE radiation_3d_data_averaging( mode, variable )
2470 
2471
2472    USE control_parameters
2473
2474    USE indices
2475
2476    USE kinds
2477
2478    IMPLICIT NONE
2479
2480    CHARACTER (LEN=*) ::  mode    !<
2481    CHARACTER (LEN=*) :: variable !<
2482
2483    INTEGER(iwp) ::  i !<
2484    INTEGER(iwp) ::  j !<
2485    INTEGER(iwp) ::  k !<
2486
2487    IF ( mode == 'allocate' )  THEN
2488
2489       SELECT CASE ( TRIM( variable ) )
2490
2491             CASE ( 'rad_net*' )
2492                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
2493                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
2494                ENDIF
2495                rad_net_av = 0.0_wp
2496
2497             CASE ( 'rad_lw_in' )
2498                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
2499                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2500                ENDIF
2501                rad_lw_in_av = 0.0_wp
2502
2503             CASE ( 'rad_lw_out' )
2504                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
2505                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2506                ENDIF
2507                rad_lw_out_av = 0.0_wp
2508
2509             CASE ( 'rad_lw_cs_hr' )
2510                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
2511                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2512                ENDIF
2513                rad_lw_cs_hr_av = 0.0_wp
2514
2515             CASE ( 'rad_lw_hr' )
2516                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
2517                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2518                ENDIF
2519                rad_lw_hr_av = 0.0_wp
2520
2521             CASE ( 'rad_sw_in' )
2522                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
2523                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2524                ENDIF
2525                rad_sw_in_av = 0.0_wp
2526
2527             CASE ( 'rad_sw_out' )
2528                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
2529                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2530                ENDIF
2531                rad_sw_out_av = 0.0_wp
2532
2533             CASE ( 'rad_sw_cs_hr' )
2534                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
2535                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2536                ENDIF
2537                rad_sw_cs_hr_av = 0.0_wp
2538
2539             CASE ( 'rad_sw_hr' )
2540                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
2541                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2542                ENDIF
2543                rad_sw_hr_av = 0.0_wp
2544
2545          CASE DEFAULT
2546             CONTINUE
2547
2548       END SELECT
2549
2550    ELSEIF ( mode == 'sum' )  THEN
2551
2552       SELECT CASE ( TRIM( variable ) )
2553
2554          CASE ( 'rad_net*' )
2555             DO  i = nxlg, nxrg
2556                DO  j = nysg, nyng
2557                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
2558                ENDDO
2559             ENDDO
2560
2561          CASE ( 'rad_lw_in' )
2562             DO  i = nxlg, nxrg
2563                DO  j = nysg, nyng
2564                   DO  k = nzb, nzt+1
2565                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
2566                   ENDDO
2567                ENDDO
2568             ENDDO
2569
2570          CASE ( 'rad_lw_out' )
2571             DO  i = nxlg, nxrg
2572                DO  j = nysg, nyng
2573                   DO  k = nzb, nzt+1
2574                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
2575                                             + rad_lw_out(k,j,i)
2576                   ENDDO
2577                ENDDO
2578             ENDDO
2579
2580          CASE ( 'rad_lw_cs_hr' )
2581             DO  i = nxlg, nxrg
2582                DO  j = nysg, nyng
2583                   DO  k = nzb, nzt+1
2584                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
2585                                               + rad_lw_cs_hr(k,j,i)
2586                   ENDDO
2587                ENDDO
2588             ENDDO
2589
2590          CASE ( 'rad_lw_hr' )
2591             DO  i = nxlg, nxrg
2592                DO  j = nysg, nyng
2593                   DO  k = nzb, nzt+1
2594                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
2595                                            + rad_lw_hr(k,j,i)
2596                   ENDDO
2597                ENDDO
2598             ENDDO
2599
2600          CASE ( 'rad_sw_in' )
2601             DO  i = nxlg, nxrg
2602                DO  j = nysg, nyng
2603                   DO  k = nzb, nzt+1
2604                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
2605                                            + rad_sw_in(k,j,i)
2606                   ENDDO
2607                ENDDO
2608             ENDDO
2609
2610          CASE ( 'rad_sw_out' )
2611             DO  i = nxlg, nxrg
2612                DO  j = nysg, nyng
2613                   DO  k = nzb, nzt+1
2614                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
2615                                             + rad_sw_out(k,j,i)
2616                   ENDDO
2617                ENDDO
2618             ENDDO
2619
2620          CASE ( 'rad_sw_cs_hr' )
2621             DO  i = nxlg, nxrg
2622                DO  j = nysg, nyng
2623                   DO  k = nzb, nzt+1
2624                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
2625                                               + rad_sw_cs_hr(k,j,i)
2626                   ENDDO
2627                ENDDO
2628             ENDDO
2629
2630          CASE ( 'rad_sw_hr' )
2631             DO  i = nxlg, nxrg
2632                DO  j = nysg, nyng
2633                   DO  k = nzb, nzt+1
2634                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
2635                                            + rad_sw_hr(k,j,i)
2636                   ENDDO
2637                ENDDO
2638             ENDDO
2639
2640          CASE DEFAULT
2641             CONTINUE
2642
2643       END SELECT
2644
2645    ELSEIF ( mode == 'average' )  THEN
2646
2647       SELECT CASE ( TRIM( variable ) )
2648
2649         CASE ( 'rad_net*' )
2650             DO  i = nxlg, nxrg
2651                DO  j = nysg, nyng
2652                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, &
2653                                     KIND=wp )
2654                ENDDO
2655             ENDDO
2656
2657          CASE ( 'rad_lw_in' )
2658             DO  i = nxlg, nxrg
2659                DO  j = nysg, nyng
2660                   DO  k = nzb, nzt+1
2661                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i)                &
2662                                            / REAL( average_count_3d, KIND=wp )
2663                   ENDDO
2664                ENDDO
2665             ENDDO
2666
2667          CASE ( 'rad_lw_out' )
2668             DO  i = nxlg, nxrg
2669                DO  j = nysg, nyng
2670                   DO  k = nzb, nzt+1
2671                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
2672                                             / REAL( average_count_3d, KIND=wp )
2673                   ENDDO
2674                ENDDO
2675             ENDDO
2676
2677          CASE ( 'rad_lw_cs_hr' )
2678             DO  i = nxlg, nxrg
2679                DO  j = nysg, nyng
2680                   DO  k = nzb, nzt+1
2681                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
2682                                             / REAL( average_count_3d, KIND=wp )
2683                   ENDDO
2684                ENDDO
2685             ENDDO
2686
2687          CASE ( 'rad_lw_hr' )
2688             DO  i = nxlg, nxrg
2689                DO  j = nysg, nyng
2690                   DO  k = nzb, nzt+1
2691                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
2692                                            / REAL( average_count_3d, KIND=wp )
2693                   ENDDO
2694                ENDDO
2695             ENDDO
2696
2697          CASE ( 'rad_sw_in' )
2698             DO  i = nxlg, nxrg
2699                DO  j = nysg, nyng
2700                   DO  k = nzb, nzt+1
2701                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
2702                                            / REAL( average_count_3d, KIND=wp )
2703                   ENDDO
2704                ENDDO
2705             ENDDO
2706
2707          CASE ( 'rad_sw_out' )
2708             DO  i = nxlg, nxrg
2709                DO  j = nysg, nyng
2710                   DO  k = nzb, nzt+1
2711                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
2712                                             / REAL( average_count_3d, KIND=wp )
2713                   ENDDO
2714                ENDDO
2715             ENDDO
2716
2717          CASE ( 'rad_sw_cs_hr' )
2718             DO  i = nxlg, nxrg
2719                DO  j = nysg, nyng
2720                   DO  k = nzb, nzt+1
2721                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
2722                                             / REAL( average_count_3d, KIND=wp )
2723                   ENDDO
2724                ENDDO
2725             ENDDO
2726
2727          CASE ( 'rad_sw_hr' )
2728             DO  i = nxlg, nxrg
2729                DO  j = nysg, nyng
2730                   DO  k = nzb, nzt+1
2731                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
2732                                            / REAL( average_count_3d, KIND=wp )
2733                   ENDDO
2734                ENDDO
2735             ENDDO
2736
2737       END SELECT
2738
2739    ENDIF
2740
2741END SUBROUTINE radiation_3d_data_averaging
2742
2743
2744!------------------------------------------------------------------------------!
2745!
2746! Description:
2747! ------------
2748!> Subroutine defining appropriate grid for netcdf variables.
2749!> It is called out from subroutine netcdf.
2750!------------------------------------------------------------------------------!
2751SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2752   
2753    IMPLICIT NONE
2754
2755    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
2756    LOGICAL, INTENT(OUT)           ::  found       !<
2757    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
2758    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
2759    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
2760
2761    found  = .TRUE.
2762
2763
2764!
2765!-- Check for the grid
2766    SELECT CASE ( TRIM( var ) )
2767
2768       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
2769              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
2770              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
2771              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
2772              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
2773          grid_x = 'x'
2774          grid_y = 'y'
2775          grid_z = 'zu'
2776
2777       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
2778              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
2779              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
2780              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
2781          grid_x = 'x'
2782          grid_y = 'y'
2783          grid_z = 'zw'
2784
2785
2786       CASE DEFAULT
2787          found  = .FALSE.
2788          grid_x = 'none'
2789          grid_y = 'none'
2790          grid_z = 'none'
2791
2792        END SELECT
2793
2794    END SUBROUTINE radiation_define_netcdf_grid
2795
2796!------------------------------------------------------------------------------!
2797!
2798! Description:
2799! ------------
2800!> Subroutine defining 3D output variables
2801!------------------------------------------------------------------------------!
2802 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
2803                                      local_pf, two_d )
2804 
2805    USE indices
2806
2807    USE kinds
2808
2809
2810    IMPLICIT NONE
2811
2812    CHARACTER (LEN=*) ::  grid     !<
2813    CHARACTER (LEN=*) ::  mode     !<
2814    CHARACTER (LEN=*) ::  variable !<
2815
2816    INTEGER(iwp) ::  av !<
2817    INTEGER(iwp) ::  i  !<
2818    INTEGER(iwp) ::  j  !<
2819    INTEGER(iwp) ::  k  !<
2820
2821    LOGICAL      ::  found !<
2822    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2823
2824    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
2825
2826    found = .TRUE.
2827
2828    SELECT CASE ( TRIM( variable ) )
2829
2830       CASE ( 'rad_net*_xy' )        ! 2d-array
2831          IF ( av == 0 ) THEN
2832             DO  i = nxl, nxr
2833                DO  j = nys, nyn
2834                   local_pf(i,j,nzb+1) = rad_net(j,i)
2835                ENDDO
2836             ENDDO
2837          ELSE
2838             DO  i = nxl, nxr
2839                DO  j = nys, nyn 
2840                   local_pf(i,j,nzb+1) = rad_net_av(j,i)
2841                ENDDO
2842             ENDDO
2843          ENDIF
2844          two_d = .TRUE.
2845          grid = 'zu1'
2846
2847 
2848       CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
2849          IF ( av == 0 ) THEN
2850             DO  i = nxl, nxr
2851                DO  j = nys, nyn
2852                   DO  k = nzb, nzt+1
2853                      local_pf(i,j,k) = rad_lw_in(k,j,i)
2854                   ENDDO
2855                ENDDO
2856             ENDDO
2857          ELSE
2858             DO  i = nxl, nxr
2859                DO  j = nys, nyn 
2860                   DO  k = nzb, nzt+1
2861                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
2862                   ENDDO
2863                ENDDO
2864             ENDDO
2865          ENDIF
2866          IF ( mode == 'xy' )  grid = 'zu'
2867
2868       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
2869          IF ( av == 0 ) THEN
2870             DO  i = nxl, nxr
2871                DO  j = nys, nyn
2872                   DO  k = nzb, nzt+1
2873                      local_pf(i,j,k) = rad_lw_out(k,j,i)
2874                   ENDDO
2875                ENDDO
2876             ENDDO
2877          ELSE
2878             DO  i = nxl, nxr
2879                DO  j = nys, nyn 
2880                   DO  k = nzb, nzt+1
2881                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
2882                   ENDDO
2883                ENDDO
2884             ENDDO
2885          ENDIF   
2886          IF ( mode == 'xy' )  grid = 'zu'
2887
2888       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
2889          IF ( av == 0 ) THEN
2890             DO  i = nxl, nxr
2891                DO  j = nys, nyn
2892                   DO  k = nzb, nzt+1
2893                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
2894                   ENDDO
2895                ENDDO
2896             ENDDO
2897          ELSE
2898             DO  i = nxl, nxr
2899                DO  j = nys, nyn 
2900                   DO  k = nzb, nzt+1
2901                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
2902                   ENDDO
2903                ENDDO
2904             ENDDO
2905          ENDIF
2906          IF ( mode == 'xy' )  grid = 'zw'
2907
2908       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
2909          IF ( av == 0 ) THEN
2910             DO  i = nxl, nxr
2911                DO  j = nys, nyn
2912                   DO  k = nzb, nzt+1
2913                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
2914                   ENDDO
2915                ENDDO
2916             ENDDO
2917          ELSE
2918             DO  i = nxl, nxr
2919                DO  j = nys, nyn 
2920                   DO  k = nzb, nzt+1
2921                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
2922                   ENDDO
2923                ENDDO
2924             ENDDO
2925          ENDIF
2926          IF ( mode == 'xy' )  grid = 'zw'
2927
2928       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
2929          IF ( av == 0 ) THEN
2930             DO  i = nxl, nxr
2931                DO  j = nys, nyn
2932                   DO  k = nzb, nzt+1
2933                      local_pf(i,j,k) = rad_sw_in(k,j,i)
2934                   ENDDO
2935                ENDDO
2936             ENDDO
2937          ELSE
2938             DO  i = nxl, nxr
2939                DO  j = nys, nyn 
2940                   DO  k = nzb, nzt+1
2941                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2942                   ENDDO
2943                ENDDO
2944             ENDDO
2945          ENDIF
2946          IF ( mode == 'xy' )  grid = 'zu'
2947
2948       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
2949          IF ( av == 0 ) THEN
2950             DO  i = nxl, nxr
2951                DO  j = nys, nyn
2952                   DO  k = nzb, nzt+1
2953                      local_pf(i,j,k) = rad_sw_out(k,j,i)
2954                   ENDDO
2955                ENDDO
2956             ENDDO
2957          ELSE
2958             DO  i = nxl, nxr
2959                DO  j = nys, nyn 
2960                   DO  k = nzb, nzt+1
2961                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2962                   ENDDO
2963                ENDDO
2964             ENDDO
2965          ENDIF
2966          IF ( mode == 'xy' )  grid = 'zu'
2967
2968       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
2969          IF ( av == 0 ) THEN
2970             DO  i = nxl, nxr
2971                DO  j = nys, nyn
2972                   DO  k = nzb, nzt+1
2973                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2974                   ENDDO
2975                ENDDO
2976             ENDDO
2977          ELSE
2978             DO  i = nxl, nxr
2979                DO  j = nys, nyn 
2980                   DO  k = nzb, nzt+1
2981                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2982                   ENDDO
2983                ENDDO
2984             ENDDO
2985          ENDIF
2986          IF ( mode == 'xy' )  grid = 'zw'
2987
2988       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
2989          IF ( av == 0 ) THEN
2990             DO  i = nxl, nxr
2991                DO  j = nys, nyn
2992                   DO  k = nzb, nzt+1
2993                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
2994                   ENDDO
2995                ENDDO
2996             ENDDO
2997          ELSE
2998             DO  i = nxl, nxr
2999                DO  j = nys, nyn 
3000                   DO  k = nzb, nzt+1
3001                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3002                   ENDDO
3003                ENDDO
3004             ENDDO
3005          ENDIF
3006          IF ( mode == 'xy' )  grid = 'zw'
3007
3008       CASE DEFAULT
3009          found = .FALSE.
3010          grid  = 'none'
3011
3012    END SELECT
3013 
3014 END SUBROUTINE radiation_data_output_2d
3015
3016
3017!------------------------------------------------------------------------------!
3018!
3019! Description:
3020! ------------
3021!> Subroutine defining 3D output variables
3022!------------------------------------------------------------------------------!
3023 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
3024 
3025
3026    USE indices
3027
3028    USE kinds
3029
3030
3031    IMPLICIT NONE
3032
3033    CHARACTER (LEN=*) ::  variable !<
3034
3035    INTEGER(iwp) ::  av    !<
3036    INTEGER(iwp) ::  i     !<
3037    INTEGER(iwp) ::  j     !<
3038    INTEGER(iwp) ::  k     !<
3039
3040    LOGICAL      ::  found !<
3041
3042    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
3043
3044
3045    found = .TRUE.
3046
3047
3048    SELECT CASE ( TRIM( variable ) )
3049
3050      CASE ( 'rad_sw_in' )
3051         IF ( av == 0 )  THEN
3052            DO  i = nxl, nxr
3053               DO  j = nys, nyn
3054                  DO  k = nzb, nzt+1
3055                     local_pf(i,j,k) = rad_sw_in(k,j,i)
3056                  ENDDO
3057               ENDDO
3058            ENDDO
3059         ELSE
3060            DO  i = nxl, nxr
3061               DO  j = nys, nyn
3062                  DO  k = nzb, nzt+1
3063                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
3064                  ENDDO
3065               ENDDO
3066            ENDDO
3067         ENDIF
3068
3069      CASE ( 'rad_sw_out' )
3070         IF ( av == 0 )  THEN
3071            DO  i = nxl, nxr
3072               DO  j = nys, nyn
3073                  DO  k = nzb, nzt+1
3074                     local_pf(i,j,k) = rad_sw_out(k,j,i)
3075                  ENDDO
3076               ENDDO
3077            ENDDO
3078         ELSE
3079            DO  i = nxl, nxr
3080               DO  j = nys, nyn
3081                  DO  k = nzb, nzt+1
3082                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
3083                  ENDDO
3084               ENDDO
3085            ENDDO
3086         ENDIF
3087
3088      CASE ( 'rad_sw_cs_hr' )
3089         IF ( av == 0 )  THEN
3090            DO  i = nxl, nxr
3091               DO  j = nys, nyn
3092                  DO  k = nzb, nzt+1
3093                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
3094                  ENDDO
3095               ENDDO
3096            ENDDO
3097         ELSE
3098            DO  i = nxl, nxr
3099               DO  j = nys, nyn
3100                  DO  k = nzb, nzt+1
3101                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
3102                  ENDDO
3103               ENDDO
3104            ENDDO
3105         ENDIF
3106
3107      CASE ( 'rad_sw_hr' )
3108         IF ( av == 0 )  THEN
3109            DO  i = nxl, nxr
3110               DO  j = nys, nyn
3111                  DO  k = nzb, nzt+1
3112                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
3113                  ENDDO
3114               ENDDO
3115            ENDDO
3116         ELSE
3117            DO  i = nxl, nxr
3118               DO  j = nys, nyn
3119                  DO  k = nzb, nzt+1
3120                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3121                  ENDDO
3122               ENDDO
3123            ENDDO
3124         ENDIF
3125
3126      CASE ( 'rad_lw_in' )
3127         IF ( av == 0 )  THEN
3128            DO  i = nxl, nxr
3129               DO  j = nys, nyn
3130                  DO  k = nzb, nzt+1
3131                     local_pf(i,j,k) = rad_lw_in(k,j,i)
3132                  ENDDO
3133               ENDDO
3134            ENDDO
3135         ELSE
3136            DO  i = nxl, nxr
3137               DO  j = nys, nyn
3138                  DO  k = nzb, nzt+1
3139                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
3140                  ENDDO
3141               ENDDO
3142            ENDDO
3143         ENDIF
3144
3145      CASE ( 'rad_lw_out' )
3146         IF ( av == 0 )  THEN
3147            DO  i = nxl, nxr
3148               DO  j = nys, nyn
3149                  DO  k = nzb, nzt+1
3150                     local_pf(i,j,k) = rad_lw_out(k,j,i)
3151                  ENDDO
3152               ENDDO
3153            ENDDO
3154         ELSE
3155            DO  i = nxl, nxr
3156               DO  j = nys, nyn
3157                  DO  k = nzb, nzt+1
3158                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
3159                  ENDDO
3160               ENDDO
3161            ENDDO
3162         ENDIF
3163
3164      CASE ( 'rad_lw_cs_hr' )
3165         IF ( av == 0 )  THEN
3166            DO  i = nxl, nxr
3167               DO  j = nys, nyn
3168                  DO  k = nzb, nzt+1
3169                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
3170                  ENDDO
3171               ENDDO
3172            ENDDO
3173         ELSE
3174            DO  i = nxl, nxr
3175               DO  j = nys, nyn
3176                  DO  k = nzb, nzt+1
3177                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
3178                  ENDDO
3179               ENDDO
3180            ENDDO
3181         ENDIF
3182
3183      CASE ( 'rad_lw_hr' )
3184         IF ( av == 0 )  THEN
3185            DO  i = nxl, nxr
3186               DO  j = nys, nyn
3187                  DO  k = nzb, nzt+1
3188                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
3189                  ENDDO
3190               ENDDO
3191            ENDDO
3192         ELSE
3193            DO  i = nxl, nxr
3194               DO  j = nys, nyn
3195                  DO  k = nzb, nzt+1
3196                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
3197                  ENDDO
3198               ENDDO
3199            ENDDO
3200         ENDIF
3201
3202       CASE DEFAULT
3203          found = .FALSE.
3204
3205    END SELECT
3206
3207
3208 END SUBROUTINE radiation_data_output_3d
3209
3210!------------------------------------------------------------------------------!
3211!
3212! Description:
3213! ------------
3214!> Subroutine defining masked data output
3215!------------------------------------------------------------------------------!
3216 SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
3217 
3218    USE control_parameters
3219       
3220    USE indices
3221   
3222    USE kinds
3223   
3224
3225    IMPLICIT NONE
3226
3227    CHARACTER (LEN=*) ::  variable   !<
3228
3229    INTEGER(iwp) ::  av   !<
3230    INTEGER(iwp) ::  i    !<
3231    INTEGER(iwp) ::  j    !<
3232    INTEGER(iwp) ::  k    !<
3233
3234    LOGICAL ::  found     !<
3235
3236    REAL(wp),                                                                  &
3237       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
3238          local_pf   !<
3239
3240
3241    found = .TRUE.
3242
3243    SELECT CASE ( TRIM( variable ) )
3244
3245
3246       CASE ( 'rad_lw_in' )
3247          IF ( av == 0 )  THEN
3248             DO  i = 1, mask_size_l(mid,1)
3249                DO  j = 1, mask_size_l(mid,2)
3250                   DO  k = 1, mask_size_l(mid,3)
3251                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
3252                                            mask_j(mid,j),mask_i(mid,i))
3253                    ENDDO
3254                 ENDDO
3255              ENDDO
3256          ELSE
3257             DO  i = 1, mask_size_l(mid,1)
3258                DO  j = 1, mask_size_l(mid,2)
3259                   DO  k = 1, mask_size_l(mid,3)
3260                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
3261                                               mask_j(mid,j),mask_i(mid,i))
3262                   ENDDO
3263                ENDDO
3264             ENDDO
3265          ENDIF
3266
3267       CASE ( 'rad_lw_out' )
3268          IF ( av == 0 )  THEN
3269             DO  i = 1, mask_size_l(mid,1)
3270                DO  j = 1, mask_size_l(mid,2)
3271                   DO  k = 1, mask_size_l(mid,3)
3272                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
3273                                            mask_j(mid,j),mask_i(mid,i))
3274                    ENDDO
3275                 ENDDO
3276              ENDDO
3277          ELSE
3278             DO  i = 1, mask_size_l(mid,1)
3279                DO  j = 1, mask_size_l(mid,2)
3280                   DO  k = 1, mask_size_l(mid,3)
3281                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
3282                                               mask_j(mid,j),mask_i(mid,i))
3283                   ENDDO
3284                ENDDO
3285             ENDDO
3286          ENDIF
3287
3288       CASE ( 'rad_lw_cs_hr' )
3289          IF ( av == 0 )  THEN
3290             DO  i = 1, mask_size_l(mid,1)
3291                DO  j = 1, mask_size_l(mid,2)
3292                   DO  k = 1, mask_size_l(mid,3)
3293                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
3294                                            mask_j(mid,j),mask_i(mid,i))
3295                    ENDDO
3296                 ENDDO
3297              ENDDO
3298          ELSE
3299             DO  i = 1, mask_size_l(mid,1)
3300                DO  j = 1, mask_size_l(mid,2)
3301                   DO  k = 1, mask_size_l(mid,3)
3302                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
3303                                               mask_j(mid,j),mask_i(mid,i))
3304                   ENDDO
3305                ENDDO
3306             ENDDO
3307          ENDIF
3308
3309       CASE ( 'rad_lw_hr' )
3310          IF ( av == 0 )  THEN
3311             DO  i = 1, mask_size_l(mid,1)
3312                DO  j = 1, mask_size_l(mid,2)
3313                   DO  k = 1, mask_size_l(mid,3)
3314                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
3315                                            mask_j(mid,j),mask_i(mid,i))
3316                    ENDDO
3317                 ENDDO
3318              ENDDO
3319          ELSE
3320             DO  i = 1, mask_size_l(mid,1)
3321                DO  j = 1, mask_size_l(mid,2)
3322                   DO  k = 1, mask_size_l(mid,3)
3323                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
3324                                               mask_j(mid,j),mask_i(mid,i))
3325                   ENDDO
3326                ENDDO
3327             ENDDO
3328          ENDIF
3329
3330       CASE ( 'rad_sw_in' )
3331          IF ( av == 0 )  THEN
3332             DO  i = 1, mask_size_l(mid,1)
3333                DO  j = 1, mask_size_l(mid,2)
3334                   DO  k = 1, mask_size_l(mid,3)
3335                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
3336                                            mask_j(mid,j),mask_i(mid,i))
3337                    ENDDO
3338                 ENDDO
3339              ENDDO
3340          ELSE
3341             DO  i = 1, mask_size_l(mid,1)
3342                DO  j = 1, mask_size_l(mid,2)
3343                   DO  k = 1, mask_size_l(mid,3)
3344                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
3345                                               mask_j(mid,j),mask_i(mid,i))
3346                   ENDDO
3347                ENDDO
3348             ENDDO
3349          ENDIF
3350
3351       CASE ( 'rad_sw_out' )
3352          IF ( av == 0 )  THEN
3353             DO  i = 1, mask_size_l(mid,1)
3354                DO  j = 1, mask_size_l(mid,2)
3355                   DO  k = 1, mask_size_l(mid,3)
3356                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
3357                                            mask_j(mid,j),mask_i(mid,i))
3358                    ENDDO
3359                 ENDDO
3360              ENDDO
3361          ELSE
3362             DO  i = 1, mask_size_l(mid,1)
3363                DO  j = 1, mask_size_l(mid,2)
3364                   DO  k = 1, mask_size_l(mid,3)
3365                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
3366                                               mask_j(mid,j),mask_i(mid,i))
3367                   ENDDO
3368                ENDDO
3369             ENDDO
3370          ENDIF
3371
3372       CASE ( 'rad_sw_cs_hr' )
3373          IF ( av == 0 )  THEN
3374             DO  i = 1, mask_size_l(mid,1)
3375                DO  j = 1, mask_size_l(mid,2)
3376                   DO  k = 1, mask_size_l(mid,3)
3377                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
3378                                            mask_j(mid,j),mask_i(mid,i))
3379                    ENDDO
3380                 ENDDO
3381              ENDDO
3382          ELSE
3383             DO  i = 1, mask_size_l(mid,1)
3384                DO  j = 1, mask_size_l(mid,2)
3385                   DO  k = 1, mask_size_l(mid,3)
3386                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
3387                                               mask_j(mid,j),mask_i(mid,i))
3388                   ENDDO
3389                ENDDO
3390             ENDDO
3391          ENDIF
3392
3393       CASE ( 'rad_sw_hr' )
3394          IF ( av == 0 )  THEN
3395             DO  i = 1, mask_size_l(mid,1)
3396                DO  j = 1, mask_size_l(mid,2)
3397                   DO  k = 1, mask_size_l(mid,3)
3398                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
3399                                            mask_j(mid,j),mask_i(mid,i))
3400                    ENDDO
3401                 ENDDO
3402              ENDDO
3403          ELSE
3404             DO  i = 1, mask_size_l(mid,1)
3405                DO  j = 1, mask_size_l(mid,2)
3406                   DO  k = 1, mask_size_l(mid,3)
3407                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
3408                                               mask_j(mid,j),mask_i(mid,i))
3409                   ENDDO
3410                ENDDO
3411             ENDDO
3412          ENDIF
3413
3414       CASE DEFAULT
3415          found = .FALSE.
3416
3417    END SELECT
3418
3419
3420 END SUBROUTINE radiation_data_output_mask
3421
3422
3423!------------------------------------------------------------------------------!
3424!
3425! Description:
3426! ------------
3427!> Subroutine defines masked output variables
3428!------------------------------------------------------------------------------!
3429 SUBROUTINE radiation_last_actions
3430 
3431
3432    USE control_parameters
3433       
3434    USE kinds
3435
3436    IMPLICIT NONE
3437
3438    IF ( write_binary )  THEN
3439       IF ( ALLOCATED( rad_net ) )  THEN
3440          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
3441       ENDIF
3442       IF ( ALLOCATED( rad_net_av ) )  THEN
3443          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
3444       ENDIF 
3445       IF ( ALLOCATED( rad_lw_in ) )  THEN
3446          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
3447       ENDIF
3448       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
3449          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
3450       ENDIF
3451       IF ( ALLOCATED( rad_lw_out ) )  THEN
3452          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out 
3453       ENDIF
3454       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
3455          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
3456       ENDIF
3457       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
3458          WRITE ( 14 )  'rad_lw_out_change_0 '
3459          WRITE ( 14 )  rad_lw_out_change_0
3460       ENDIF
3461       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
3462          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
3463       ENDIF
3464       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3465          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
3466       ENDIF
3467       IF ( ALLOCATED( rad_lw_hr ) )  THEN
3468          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
3469       ENDIF
3470       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
3471          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
3472       ENDIF
3473       IF ( ALLOCATED( rad_sw_in ) )  THEN
3474          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
3475       ENDIF
3476       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
3477          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
3478       ENDIF
3479       IF ( ALLOCATED( rad_sw_out ) )  THEN
3480          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
3481       ENDIF
3482       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
3483          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
3484       ENDIF
3485       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
3486          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
3487       ENDIF
3488       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3489          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
3490       ENDIF
3491       IF ( ALLOCATED( rad_sw_hr ) )  THEN
3492          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
3493       ENDIF
3494       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
3495          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
3496       ENDIF
3497
3498       WRITE ( 14 )  '*** end rad ***     '
3499
3500    ENDIF
3501
3502 END SUBROUTINE radiation_last_actions
3503
3504
3505SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,          &
3506                                        nxr_on_file, nynfa, nyn_on_file, nysfa,&
3507                                        nys_on_file, offset_xa, offset_ya,     &
3508                                        overlap_count, tmp_2d, tmp_3d )
3509 
3510
3511    USE control_parameters
3512       
3513    USE indices
3514   
3515    USE kinds
3516   
3517    USE pegrid
3518
3519    IMPLICIT NONE
3520
3521    CHARACTER (LEN=20) :: field_char   !<
3522
3523    INTEGER(iwp) ::  i               !<
3524    INTEGER(iwp) ::  k               !<
3525    INTEGER(iwp) ::  nxlc            !<
3526    INTEGER(iwp) ::  nxlf            !<
3527    INTEGER(iwp) ::  nxl_on_file     !<
3528    INTEGER(iwp) ::  nxrc            !<
3529    INTEGER(iwp) ::  nxrf            !<
3530    INTEGER(iwp) ::  nxr_on_file     !<
3531    INTEGER(iwp) ::  nync            !<
3532    INTEGER(iwp) ::  nynf            !<
3533    INTEGER(iwp) ::  nyn_on_file     !<
3534    INTEGER(iwp) ::  nysc            !<
3535    INTEGER(iwp) ::  nysf            !<
3536    INTEGER(iwp) ::  nys_on_file     !<
3537    INTEGER(iwp) ::  overlap_count   !<
3538
3539    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
3540    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
3541    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
3542    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
3543    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
3544    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
3545
3546    REAL(wp),                                                                  &
3547       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3548          tmp_2d   !<
3549
3550    REAL(wp),                                                                  &
3551       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3552          tmp_3d   !<
3553
3554    REAL(wp),                                                                  &
3555       DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3556          tmp_3d2   !<
3557
3558
3559
3560   IF ( initializing_actions == 'read_restart_data' )  THEN
3561      READ ( 13 )  field_char
3562
3563      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
3564
3565         DO  k = 1, overlap_count
3566
3567            nxlf = nxlfa(i,k)
3568            nxlc = nxlfa(i,k) + offset_xa(i,k)
3569            nxrf = nxrfa(i,k)
3570            nxrc = nxrfa(i,k) + offset_xa(i,k)
3571            nysf = nysfa(i,k)
3572            nysc = nysfa(i,k) + offset_ya(i,k)
3573            nynf = nynfa(i,k)
3574            nync = nynfa(i,k) + offset_ya(i,k)
3575
3576
3577            SELECT CASE ( TRIM( field_char ) )
3578
3579                CASE ( 'rad_net' )
3580                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
3581                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
3582                   ENDIF 
3583                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3584                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =         &
3585                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3586
3587                CASE ( 'rad_net_av' )
3588                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
3589                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
3590                   ENDIF 
3591                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3592                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
3593                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3594                CASE ( 'rad_lw_in' )
3595                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
3596                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3597                           radiation_scheme == 'constant')  THEN
3598                         ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
3599                      ELSE
3600                         ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3601                      ENDIF
3602                   ENDIF 
3603                   IF ( k == 1 )  THEN
3604                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3605                           radiation_scheme == 'constant')  THEN
3606                         READ ( 13 )  tmp_3d2
3607                         rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3608                            tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3609                      ELSE
3610                         READ ( 13 )  tmp_3d
3611                         rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3612                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3613                      ENDIF
3614                   ENDIF
3615
3616                CASE ( 'rad_lw_in_av' )
3617                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
3618                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3619                           radiation_scheme == 'constant')  THEN
3620                         ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3621                      ELSE
3622                         ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3623                      ENDIF
3624                   ENDIF 
3625                   IF ( k == 1 )  THEN
3626                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3627                           radiation_scheme == 'constant')  THEN
3628                         READ ( 13 )  tmp_3d2
3629                         rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3630                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3631                      ELSE
3632                         READ ( 13 )  tmp_3d
3633                         rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3634                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3635                      ENDIF
3636                   ENDIF
3637
3638                CASE ( 'rad_lw_out' )
3639                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
3640                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3641                           radiation_scheme == 'constant')  THEN
3642                         ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
3643                      ELSE
3644                         ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3645                      ENDIF
3646                   ENDIF 
3647                   IF ( k == 1 )  THEN
3648                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3649                           radiation_scheme == 'constant')  THEN
3650                         READ ( 13 )  tmp_3d2
3651                         rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3652                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3653                      ELSE
3654                         READ ( 13 )  tmp_3d
3655                         rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3656                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3657                      ENDIF
3658                   ENDIF
3659
3660                CASE ( 'rad_lw_out_av' )
3661                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
3662                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3663                           radiation_scheme == 'constant')  THEN
3664                         ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3665                      ELSE
3666                         ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3667                      ENDIF
3668                   ENDIF 
3669                   IF ( k == 1 )  THEN
3670                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3671                           radiation_scheme == 'constant')  THEN
3672                         READ ( 13 )  tmp_3d2
3673                         rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3674                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3675                      ELSE
3676                         READ ( 13 )  tmp_3d
3677                         rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3678                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3679                      ENDIF
3680                   ENDIF
3681
3682                CASE ( 'rad_lw_out_change_0' )
3683                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
3684                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
3685                   ENDIF
3686                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3687                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
3688                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3689
3690                CASE ( 'rad_lw_cs_hr' )
3691                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
3692                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3693                   ENDIF
3694                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3695                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3696                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3697
3698                CASE ( 'rad_lw_cs_hr_av' )
3699                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3700                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3701                   ENDIF
3702                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3703                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3704                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3705
3706                CASE ( 'rad_lw_hr' )
3707                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
3708                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3709                   ENDIF
3710                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3711                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3712                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3713
3714                CASE ( 'rad_lw_hr_av' )
3715                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
3716                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3717                   ENDIF
3718                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3719                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3720                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3721
3722                CASE ( 'rad_sw_in' )
3723                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
3724                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3725                           radiation_scheme == 'constant')  THEN
3726                         ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
3727                      ELSE
3728                         ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3729                      ENDIF
3730                   ENDIF 
3731                   IF ( k == 1 )  THEN
3732                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3733                           radiation_scheme == 'constant')  THEN
3734                         READ ( 13 )  tmp_3d2
3735                         rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3736                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3737                      ELSE
3738                         READ ( 13 )  tmp_3d
3739                         rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3740                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3741                      ENDIF
3742                   ENDIF
3743
3744                CASE ( 'rad_sw_in_av' )
3745                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
3746                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3747                           radiation_scheme == 'constant')  THEN
3748                         ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3749                      ELSE
3750                         ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3751                      ENDIF
3752                   ENDIF 
3753                   IF ( k == 1 )  THEN
3754                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3755                           radiation_scheme == 'constant')  THEN
3756                         READ ( 13 )  tmp_3d2
3757                         rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3758                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3759                      ELSE
3760                         READ ( 13 )  tmp_3d
3761                         rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3762                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3763                      ENDIF
3764                   ENDIF
3765
3766                CASE ( 'rad_sw_out' )
3767                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
3768                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3769                           radiation_scheme == 'constant')  THEN
3770                         ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
3771                      ELSE
3772                         ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3773                      ENDIF
3774                   ENDIF 
3775                   IF ( k == 1 )  THEN
3776                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3777                           radiation_scheme == 'constant')  THEN
3778                         READ ( 13 )  tmp_3d2
3779                         rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3780                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3781                      ELSE
3782                         READ ( 13 )  tmp_3d
3783                         rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3784                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3785                      ENDIF
3786                   ENDIF
3787
3788                CASE ( 'rad_sw_out_av' )
3789                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
3790                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3791                           radiation_scheme == 'constant')  THEN
3792                         ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3793                      ELSE
3794                         ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3795                      ENDIF
3796                   ENDIF 
3797                   IF ( k == 1 )  THEN
3798                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3799                           radiation_scheme == 'constant')  THEN
3800                         READ ( 13 )  tmp_3d2
3801                         rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3802                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3803                      ELSE
3804                         READ ( 13 )  tmp_3d
3805                         rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3806                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3807                      ENDIF
3808                   ENDIF
3809
3810                CASE ( 'rad_sw_cs_hr' )
3811                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
3812                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3813                   ENDIF
3814                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3815                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3816                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3817
3818                CASE ( 'rad_sw_cs_hr_av' )
3819                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3820                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3821                   ENDIF
3822                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3823                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3824                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3825
3826                CASE ( 'rad_sw_hr' )
3827                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
3828                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3829                   ENDIF
3830                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3831                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3832                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3833
3834                CASE ( 'rad_sw_hr_av' )
3835                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
3836                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3837                   ENDIF
3838                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3839                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3840                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3841
3842               CASE DEFAULT
3843                  WRITE( message_string, * ) 'unknown variable named "',       &
3844                                        TRIM( field_char ), '" found in',      &
3845                                        '&data from prior run on PE ', myid
3846                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, &
3847                                0, 6, 0 )
3848
3849            END SELECT
3850
3851         ENDDO
3852
3853         READ ( 13 )  field_char
3854
3855      ENDDO
3856   ENDIF
3857
3858 END SUBROUTINE radiation_read_restart_data
3859
3860
3861 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.