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

Last change on this file since 2011 was 2011, checked in by kanani, 5 years ago

changes related to steering and formating of urban surface model

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