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

Last change on this file since 2299 was 2299, checked in by maronga, 5 years ago

improvements for spinup mechanism

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