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

Last change on this file since 2258 was 2249, checked in by sward, 7 years ago

last commit documented

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