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

Last change on this file since 2001 was 2001, checked in by knoop, 5 years ago

last commit documented

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