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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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