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

Last change on this file since 2320 was 2318, checked in by suehring, 7 years ago

get topograpyhy top index via function call

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