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

Last change on this file since 2201 was 2201, checked in by suehring, 6 years ago

last commit documented

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