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

Last change on this file since 2328 was 2328, checked in by maronga, 7 years ago

some improvements in land surface model

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