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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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