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

Last change on this file since 1849 was 1849, checked in by hoffmann, 5 years ago

lpm_droplet_condensation improved, microphysics partially modularized

  • Property svn:keywords set to Id
File size: 87.6 KB
Line 
1!> @file radiation_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21! Adapted for modularization of microphysics
22!
23! Former revisions:
24! -----------------
25! $Id: radiation_model_mod.f90 1849 2016-04-08 11:33:18Z hoffmann $
26!
27! 1826 2016-04-07 12:01:39Z maronga
28! Further modularization.
29!
30! 1788 2016-03-10 11:01:04Z maronga
31! Added new albedo class for pavements / roads.
32!
33! 1783 2016-03-06 18:36:17Z raasch
34! palm-netcdf-module removed in order to avoid a circular module dependency,
35! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
36! added
37!
38! 1757 2016-02-22 15:49:32Z maronga
39! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
40! profiles for pressure and temperature above the LES domain.
41!
42! 1709 2015-11-04 14:47:01Z maronga
43! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
44! corrections
45!
46! 1701 2015-11-02 07:43:04Z maronga
47! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
48!
49! 1691 2015-10-26 16:17:44Z maronga
50! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
51! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
52! Added output of radiative heating rates.
53!
54! 1682 2015-10-07 23:56:08Z knoop
55! Code annotations made doxygen readable
56!
57! 1606 2015-06-29 10:43:37Z maronga
58! Added preprocessor directive __netcdf to allow for compiling without netCDF.
59! Note, however, that RRTMG cannot be used without netCDF.
60!
61! 1590 2015-05-08 13:56:27Z maronga
62! Bugfix: definition of character strings requires same length for all elements
63!
64! 1587 2015-05-04 14:19:01Z maronga
65! Added albedo class for snow
66!
67! 1585 2015-04-30 07:05:52Z maronga
68! Added support for RRTMG
69!
70! 1571 2015-03-12 16:12:49Z maronga
71! Added missing KIND attribute. Removed upper-case variable names
72!
73! 1551 2015-03-03 14:18:16Z maronga
74! Added support for data output. Various variables have been renamed. Added
75! interface for different radiation schemes (currently: clear-sky, constant, and
76! RRTM (not yet implemented).
77!
78! 1496 2014-12-02 17:25:50Z maronga
79! Initial revision
80!
81!
82! Description:
83! ------------
84!> Radiation models and interfaces
85!> @todo move variable definitions used in radiation_init only to the subroutine
86!>       as they are no longer required after initialization.
87!> @todo Output of full column vertical profiles used in RRTMG
88!> @todo Output of other rrtm arrays (such as volume mixing ratios)
89!> @todo Adapt for use with topography
90!>
91!> @note Many variables have a leading dummy dimension (0:0) in order to
92!>       match the assume-size shape expected by the RRTMG model.
93!------------------------------------------------------------------------------!
94 MODULE radiation_model_mod
95 
96    USE arrays_3d,                                                             &
97        ONLY:  dzw, hyp, pt, q, ql, zw
98
99    USE cloud_parameters,                                                      &
100        ONLY:  cp, l_d_cp, rho_l
101
102    USE constants,                                                             &
103        ONLY:  pi
104
105    USE control_parameters,                                                    &
106        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
107               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
108               surface_pressure, time_since_reference_point
109
110    USE indices,                                                               &
111        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
112
113    USE kinds
114
115    USE microphysics_mod,                                                      &
116        ONLY:  nc_const, sigma_gc
117
118#if defined ( __netcdf )
119    USE NETCDF
120#endif
121
122#if defined ( __rrtmg )
123
124    USE parrrsw,                                                               &
125        ONLY:  naerec, nbndsw
126
127    USE parrrtm,                                                               &
128        ONLY:  nbndlw
129
130    USE rrtmg_lw_init,                                                         &
131        ONLY:  rrtmg_lw_ini
132
133    USE rrtmg_sw_init,                                                         &
134        ONLY:  rrtmg_sw_ini
135
136    USE rrtmg_lw_rad,                                                          &
137        ONLY:  rrtmg_lw
138
139    USE rrtmg_sw_rad,                                                          &
140        ONLY:  rrtmg_sw
141#endif
142
143
144
145    IMPLICIT NONE
146
147    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
148
149!
150!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
151    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
152                                   'user defined                         ', & !  0
153                                   'ocean                                ', & !  1
154                                   'mixed farming, tall grassland        ', & !  2
155                                   'tall/medium grassland                ', & !  3
156                                   'evergreen shrubland                  ', & !  4
157                                   'short grassland/meadow/shrubland     ', & !  5
158                                   'evergreen needleleaf forest          ', & !  6
159                                   'mixed deciduous evergreen forest     ', & !  7
160                                   'deciduous forest                     ', & !  8
161                                   'tropical evergreen broadleaved forest', & !  9
162                                   'medium/tall grassland/woodland       ', & ! 10
163                                   'desert, sandy                        ', & ! 11
164                                   'desert, rocky                        ', & ! 12
165                                   'tundra                               ', & ! 13
166                                   'land ice                             ', & ! 14
167                                   'sea ice                              ', & ! 15
168                                   'snow                                 ', & ! 16
169                                   'pavement/roads                       '  & ! 17
170                                                         /)
171
172    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
173                    day,                 & !< current day of the year
174                    day_init     = 172,  & !< day of the year at model start (21/06)
175                    dots_rad     = 0       !< starting index for timeseries output
176
177    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
178                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
179                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
180                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
181                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
182                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
183                sw_radiation = .TRUE.                   !< flag parameter indicing whether shortwave radiation shall be calculated
184
185
186    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
187                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
188                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
189                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
190
191    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
192                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
193                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
194                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
195                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
196                decl_1,                          & !< declination coef. 1
197                decl_2,                          & !< declination coef. 2
198                decl_3,                          & !< declination coef. 3
199                dt_radiation = 0.0_wp,           & !< radiation model timestep
200                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
201                lambda = 0.0_wp,                 & !< longitude in degrees
202                lon = 0.0_wp,                    & !< longitude in radians
203                lat = 0.0_wp,                    & !< latitude in radians
204                net_radiation = 0.0_wp,          & !< net radiation at surface
205                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
206                sky_trans,                       & !< sky transmissivity
207                time_radiation = 0.0_wp,         & !< time since last call of radiation code
208                time_utc,                        & !< current time in UTC
209                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
210
211    REAL(wp), DIMENSION(0:0) ::  zenith        !< solar zenith angle
212
213    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
214                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
215                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
216                rad_net,                     & !< net radiation at the surface
217                rad_net_av                     !< average of rad_net
218
219!
220!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
221!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
222    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
223                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
224                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
225                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
226                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
227                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
228                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
229                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
230                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
231                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
232                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
233                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
234                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
235                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
236                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
237                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
238                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
239                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
240                                 /), (/ 3, 17 /) )
241
242    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
243                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
244                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
245                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
246                        rad_lw_hr_av,                  & !< average of rad_sw_hr
247                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
248                        rad_lw_in_av,                  & !< average of rad_lw_in
249                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
250                        rad_lw_out_av,                 & !< average of rad_lw_out
251                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
252                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
253                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
254                        rad_sw_hr_av,                  & !< average of rad_sw_hr
255                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
256                        rad_sw_in_av,                  & !< average of rad_sw_in
257                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
258                        rad_sw_out_av                    !< average of rad_sw_out
259
260
261!
262!-- Variables and parameters used in RRTMG only
263#if defined ( __rrtmg )
264    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
265
266
267!
268!-- Flag parameters for RRTMGS (should not be changed)
269    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
270                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
271                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
272                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
273                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
274                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
275
276!
277!-- The following variables should be only changed with care, as this will
278!-- require further setting of some variables, which is currently not
279!-- implemented (aerosols, ice phase).
280    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
281                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
282                    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)
283                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
284
285    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
286
287    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
288
289    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
290
291    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
292                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
293                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
294                                           t_snd          !< actual temperature from sounding data (hPa)
295
296    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
297                                             aldir,          & !< longwave direct albedo solar angle of 60°
298                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
299                                             asdir,          & !< shortwave direct albedo solar angle of 60°
300                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
301                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
302                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
303                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
304                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
305                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
306                                             rrtm_cldfr,     & !< cloud fraction (0,1)
307                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
308                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
309                                             rrtm_emis,      & !< surface emissivity (0-1)   
310                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
311                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
312                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
313                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
314                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
315                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
316                                             rrtm_reice,     & !< cloud ice effective radius (microns)
317                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
318                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
319                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
320                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
321                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
322                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
323                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
324                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
325                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
326                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
327                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
328                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
329                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
330                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
331                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
332                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
333                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
334
335!
336!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
337    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
338                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
339                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
340                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
341                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
342                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
343                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
344                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
345                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
346                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
347                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
348                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
349                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
350                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
351                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
352                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
353                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
354                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
355
356#endif
357
358    INTERFACE radiation_check_data_output
359       MODULE PROCEDURE radiation_check_data_output
360    END INTERFACE radiation_check_data_output
361
362    INTERFACE radiation_check_data_output_pr
363       MODULE PROCEDURE radiation_check_data_output_pr
364    END INTERFACE radiation_check_data_output_pr
365 
366    INTERFACE radiation_check_parameters
367       MODULE PROCEDURE radiation_check_parameters
368    END INTERFACE radiation_check_parameters
369 
370    INTERFACE radiation_clearsky
371       MODULE PROCEDURE radiation_clearsky
372    END INTERFACE radiation_clearsky
373 
374    INTERFACE radiation_header
375       MODULE PROCEDURE radiation_header
376    END INTERFACE radiation_header
377 
378    INTERFACE radiation_init
379       MODULE PROCEDURE radiation_init
380    END INTERFACE radiation_init
381
382    INTERFACE radiation_parin
383       MODULE PROCEDURE radiation_parin
384    END INTERFACE radiation_parin
385   
386    INTERFACE radiation_rrtmg
387       MODULE PROCEDURE radiation_rrtmg
388    END INTERFACE radiation_rrtmg
389
390    INTERFACE radiation_tendency
391       MODULE PROCEDURE radiation_tendency
392       MODULE PROCEDURE radiation_tendency_ij
393    END INTERFACE radiation_tendency
394
395    SAVE
396
397    PRIVATE
398
399!
400!-- Public functions
401    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
402           radiation_check_parameters, radiation_clearsky, radiation_header,   &
403           radiation_init, radiation_parin, radiation_rrtmg, radiation_tendency
404   
405!
406!-- Public variables and constants
407    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
408           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
409           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
410           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
411           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
412           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
413           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
414
415
416#if defined ( __rrtmg )
417    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
418#endif
419
420 CONTAINS
421
422!------------------------------------------------------------------------------!
423! Description:
424! ------------
425!> Check data output for radiation model
426!------------------------------------------------------------------------------!
427    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
428 
429 
430       USE control_parameters,                                                 &
431           ONLY: data_output, message_string
432
433       IMPLICIT NONE
434
435       CHARACTER (LEN=*) ::  unit     !<
436       CHARACTER (LEN=*) ::  var      !<
437
438       INTEGER(iwp) :: i
439       INTEGER(iwp) :: ilen
440       INTEGER(iwp) :: k
441
442       SELECT CASE ( TRIM( var ) )
443
444         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
445                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
446             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
447                message_string = '"output of "' // TRIM( var ) // '" requi' // &
448                                 'res radiation = .TRUE. and ' //              &
449                                 'radiation_scheme = "rrtmg"'
450                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
451             ENDIF
452             unit = 'W/m2'     
453
454          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
455                 'rrtm_asdir*' )
456             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
457                message_string = 'illegal value for data_output: "' //         &
458                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
459                                 'cross sections are allowed for this value'
460                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
461             ENDIF
462             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
463                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
464                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
465                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
466                     TRIM( var ) == 'rrtm_asdir*'      )                       &
467                THEN
468                   message_string = 'output of "' // TRIM( var ) // '" require'&
469                                    // 's radiation = .TRUE. and radiation_sch'&
470                                    // 'eme = "rrtmg"'
471                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
472                ENDIF
473             ENDIF
474
475             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
476             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
477             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
478             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
479             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
480
481          CASE DEFAULT
482             unit = 'illegal'
483
484       END SELECT
485
486
487    END SUBROUTINE radiation_check_data_output
488
489!------------------------------------------------------------------------------!
490! Description:
491! ------------
492!> Check data output of profiles for radiation model
493!------------------------------------------------------------------------------! 
494    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
495 
496       USE arrays_3d,                                                          &
497           ONLY: zu
498
499       USE control_parameters,                                                 &
500           ONLY: data_output_pr, message_string
501
502       USE indices
503
504       USE profil_parameter
505
506       USE statistics
507
508       IMPLICIT NONE
509   
510       CHARACTER (LEN=*) ::  unit      !<
511       CHARACTER (LEN=*) ::  variable  !<
512       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
513 
514       INTEGER(iwp) ::  user_pr_index !<
515       INTEGER(iwp) ::  var_count     !<
516
517       SELECT CASE ( TRIM( variable ) )
518       
519         CASE ( 'rad_net' )
520             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
521             THEN
522                message_string = 'data_output_pr = ' //                        &
523                                 TRIM( data_output_pr(var_count) ) // ' is' // &
524                                 'not available for radiation = .FALSE. or ' //&
525                                 'radiation_scheme = "constant"'
526                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
527             ELSE
528                dopr_index(var_count) = 101
529                dopr_unit  = 'W/m2'
530                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
531                unit = dopr_unit
532             ENDIF
533
534          CASE ( 'rad_lw_in' )
535             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
536             THEN
537                message_string = 'data_output_pr = ' //                        &
538                                 TRIM( data_output_pr(var_count) ) // ' is' // &
539                                 'not available for radiation = .FALSE. or ' //&
540                                 'radiation_scheme = "constant"'
541                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
542             ELSE
543                dopr_index(var_count) = 102
544                dopr_unit  = 'W/m2'
545                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
546                unit = dopr_unit
547             ENDIF
548
549          CASE ( 'rad_lw_out' )
550             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
551             THEN
552                message_string = 'data_output_pr = ' //                        &
553                                 TRIM( data_output_pr(var_count) ) // ' is' // &
554                                 'not available for radiation = .FALSE. or ' //&
555                                 'radiation_scheme = "constant"'
556                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
557             ELSE
558                dopr_index(var_count) = 103
559                dopr_unit  = 'W/m2'
560                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
561                unit = dopr_unit   
562             ENDIF
563
564          CASE ( 'rad_sw_in' )
565             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
566             THEN
567                message_string = 'data_output_pr = ' //                        &
568                                 TRIM( data_output_pr(var_count) ) // ' is' // &
569                                 'not available for radiation = .FALSE. or ' //&
570                                 'radiation_scheme = "constant"'
571                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
572             ELSE
573                dopr_index(var_count) = 104
574                dopr_unit  = 'W/m2'
575                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
576                unit = dopr_unit
577             ENDIF
578
579          CASE ( 'rad_sw_out')
580             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
581             THEN
582                message_string = 'data_output_pr = ' //                        &
583                                 TRIM( data_output_pr(var_count) ) // ' is' // &
584                                 'not available for radiation = .FALSE. or ' //&
585                                 'radiation_scheme = "constant"'
586                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
587             ELSE
588                dopr_index(var_count) = 105
589                dopr_unit  = 'W/m2'
590                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
591                unit = dopr_unit
592             ENDIF
593
594          CASE ( 'rad_lw_cs_hr' )
595             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
596             THEN
597                message_string = 'data_output_pr = ' //                        &
598                                 TRIM( data_output_pr(var_count) ) // ' is' // &
599                                 'not available for radiation = .FALSE. or ' //&
600                                 'radiation_scheme /= "rrtmg"'
601                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
602             ELSE
603                dopr_index(var_count) = 106
604                dopr_unit  = 'K/h'
605                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
606                unit = dopr_unit
607             ENDIF
608
609          CASE ( 'rad_lw_hr' )
610             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
611             THEN
612                message_string = 'data_output_pr = ' //                        &
613                                 TRIM( data_output_pr(var_count) ) // ' is' // &
614                                 'not available for radiation = .FALSE. or ' //&
615                                 'radiation_scheme /= "rrtmg"'
616                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
617             ELSE
618                dopr_index(var_count) = 107
619                dopr_unit  = 'K/h'
620                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
621                unit = dopr_unit
622             ENDIF
623
624          CASE ( 'rad_sw_cs_hr' )
625             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
626             THEN
627                message_string = 'data_output_pr = ' //                        &
628                                 TRIM( data_output_pr(var_count) ) // ' is' // &
629                                 'not available for radiation = .FALSE. or ' //&
630                                 'radiation_scheme /= "rrtmg"'
631                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
632             ELSE
633                dopr_index(var_count) = 108
634                dopr_unit  = 'K/h'
635                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
636                unit = dopr_unit
637             ENDIF
638
639          CASE ( 'rad_sw_hr' )
640             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
641             THEN
642                message_string = 'data_output_pr = ' //                        &
643                                 TRIM( data_output_pr(var_count) ) // ' is' // &
644                                 'not available for radiation = .FALSE. or ' //&
645                                 'radiation_scheme /= "rrtmg"'
646                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
647             ELSE
648                dopr_index(var_count) = 109
649                dopr_unit  = 'K/h'
650                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
651                unit = dopr_unit
652             ENDIF
653
654
655          CASE DEFAULT
656             unit = 'illegal'
657
658       END SELECT
659
660
661    END SUBROUTINE radiation_check_data_output_pr
662 
663 
664!------------------------------------------------------------------------------!
665! Description:
666! ------------
667!> Check parameters routine for radiation model
668!------------------------------------------------------------------------------!
669    SUBROUTINE radiation_check_parameters
670
671       USE control_parameters,                                                 &
672           ONLY: message_string, topography
673                 
674   
675       IMPLICIT NONE
676
677       IF ( radiation_scheme /= 'constant'   .AND.                             &
678            radiation_scheme /= 'clear-sky'  .AND.                             &
679            radiation_scheme /= 'rrtmg' )  THEN
680          message_string = 'unknown radiation_scheme = '//                     &
681                           TRIM( radiation_scheme )
682          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
683       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
684#if ! defined ( __rrtmg )
685          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
686                           'compilation of PALM with pre-processor ' //        &
687                           'directive -D__rrtmg'
688          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
689#endif
690#if defined ( __rrtmg ) && ! defined( __netcdf )
691          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
692                           'the use of NetCDF (preprocessor directive ' //     &
693                           '-D__netcdf'
694          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
695#endif
696
697       ENDIF
698
699       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
700            radiation_scheme == 'clear-sky')  THEN
701          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
702                           'with albedo_type = 0 requires setting of albedo'// &
703                           ' /= 9999999.9'
704          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
705       ENDIF
706
707       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
708          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
709          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
710          ) ) THEN
711          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
712                           'with albedo_type = 0 requires setting of ' //      &
713                           'albedo_lw_dif /= 9999999.9' //                     &
714                           'albedo_lw_dir /= 9999999.9' //                     &
715                           'albedo_sw_dif /= 9999999.9 and' //                 &
716                           'albedo_sw_dir /= 9999999.9'
717          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
718       ENDIF
719
720       IF ( topography /= 'flat' )  THEN
721          message_string = 'radiation scheme cannot be used ' //               & 
722                           'in combination with  topography /= "flat"'
723          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
724       ENDIF
725 
726    END SUBROUTINE radiation_check_parameters
727 
728 
729!------------------------------------------------------------------------------!
730! Description:
731! ------------
732!> Initialization of the radiation model
733!------------------------------------------------------------------------------!
734    SUBROUTINE radiation_init
735   
736       IMPLICIT NONE
737
738!
739!--    Allocate array for storing the surface net radiation
740       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
741          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
742          rad_net = 0.0_wp
743       ENDIF
744
745!
746!--    Allocate array for storing the surface net radiation
747       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
748          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
749          rad_lw_out_change_0 = 0.0_wp
750       ENDIF
751
752!
753!--    Fix net radiation in case of radiation_scheme = 'constant'
754       IF ( radiation_scheme == 'constant' )  THEN
755          rad_net = net_radiation
756          radiation = .FALSE.
757!
758!--    Calculate orbital constants
759       ELSE
760          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
761          decl_2 = 2.0_wp * pi / 365.0_wp
762          decl_3 = decl_2 * 81.0_wp
763          lat    = phi * pi / 180.0_wp
764          lon    = lambda * pi / 180.0_wp
765       ENDIF
766
767
768       IF ( radiation_scheme == 'clear-sky' )  THEN
769
770          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
771
772          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
773             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
774          ENDIF
775          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
776             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
777          ENDIF
778
779          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
780             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
781          ENDIF
782          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
783             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
784          ENDIF
785
786          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
787             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
788          ENDIF
789          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
790             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
791          ENDIF
792
793          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
794             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
795          ENDIF
796          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
797             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
798          ENDIF
799
800          rad_sw_in  = 0.0_wp
801          rad_sw_out = 0.0_wp
802          rad_lw_in  = 0.0_wp
803          rad_lw_out = 0.0_wp
804
805!
806!--       Overwrite albedo if manually set in parameter file
807          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
808             albedo = albedo_pars(2,albedo_type)
809          ENDIF
810   
811          alpha = albedo
812 
813!
814!--    Initialization actions for RRTMG
815       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
816#if defined ( __rrtmg )
817!
818!--       Allocate albedos
819          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
820          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
821          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
822          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
823          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
824          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
825          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
826          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
827
828          IF ( albedo_type /= 0 )  THEN
829             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
830                albedo_lw_dif = albedo_pars(0,albedo_type)
831                albedo_lw_dir = albedo_lw_dif
832             ENDIF
833             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
834                albedo_sw_dif = albedo_pars(1,albedo_type)
835                albedo_sw_dir = albedo_sw_dif
836             ENDIF
837          ENDIF
838
839          aldif(:,:) = albedo_lw_dif
840          aldir(:,:) = albedo_lw_dir
841          asdif(:,:) = albedo_sw_dif
842          asdir(:,:) = albedo_sw_dir
843!
844!--       Calculate initial values of current (cosine of) the zenith angle and
845!--       whether the sun is up
846          CALL calc_zenith     
847!
848!--       Calculate initial surface albedo
849          IF ( .NOT. constant_albedo )  THEN
850             CALL calc_albedo
851          ELSE
852             rrtm_aldif(0,:,:) = aldif(:,:)
853             rrtm_aldir(0,:,:) = aldir(:,:)
854             rrtm_asdif(0,:,:) = asdif(:,:) 
855             rrtm_asdir(0,:,:) = asdir(:,:)   
856          ENDIF
857
858!
859!--       Allocate surface emissivity
860          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
861          rrtm_emis = emissivity
862
863!
864!--       Allocate 3d arrays of radiative fluxes and heating rates
865          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
866             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
867             rad_sw_in = 0.0_wp
868          ENDIF
869
870          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
871             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
872          ENDIF
873
874          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
875             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
876             rad_sw_out = 0.0_wp
877          ENDIF
878
879          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
880             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
881          ENDIF
882
883          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
884             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
885             rad_sw_hr = 0.0_wp
886          ENDIF
887
888          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
889             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
890             rad_sw_hr_av = 0.0_wp
891          ENDIF
892
893          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
894             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
895             rad_sw_cs_hr = 0.0_wp
896          ENDIF
897
898          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
899             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
900             rad_sw_cs_hr_av = 0.0_wp
901          ENDIF
902
903          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
904             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
905             rad_lw_in     = 0.0_wp
906          ENDIF
907
908          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
909             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
910          ENDIF
911
912          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
913             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
914            rad_lw_out    = 0.0_wp
915          ENDIF
916
917          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
918             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
919          ENDIF
920
921          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
922             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
923             rad_lw_hr = 0.0_wp
924          ENDIF
925
926          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
927             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
928             rad_lw_hr_av = 0.0_wp
929          ENDIF
930
931          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
932             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
933             rad_lw_cs_hr = 0.0_wp
934          ENDIF
935
936          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
937             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
938             rad_lw_cs_hr_av = 0.0_wp
939          ENDIF
940
941          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
942          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
943          rad_sw_cs_in  = 0.0_wp
944          rad_sw_cs_out = 0.0_wp
945
946          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
947          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
948          rad_lw_cs_in  = 0.0_wp
949          rad_lw_cs_out = 0.0_wp
950
951!
952!--       Allocate dummy array for storing surface temperature
953          ALLOCATE ( rrtm_tsfc(1) )
954
955!
956!--       Initialize RRTMG
957          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
958          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
959
960!
961!--       Set input files for RRTMG
962          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
963          IF ( .NOT. snd_exists )  THEN
964             rrtm_input_file = "rrtmg_lw.nc"
965          ENDIF
966
967!
968!--       Read vertical layers for RRTMG from sounding data
969!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
970!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
971!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
972          CALL read_sounding_data
973
974!
975!--       Read trace gas profiles from file. This routine provides
976!--       the rrtm_ arrays (1:nzt_rad+1)
977          CALL read_trace_gas_data
978#endif
979       ENDIF
980
981!
982!--    Perform user actions if required
983       CALL user_init_radiation
984
985!
986!--    Calculate radiative fluxes at model start
987       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
988          IF ( radiation_scheme == 'clear-sky' )  THEN
989             CALL radiation_clearsky
990          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
991             CALL radiation_rrtmg
992          ENDIF
993       ENDIF
994
995       RETURN
996
997    END SUBROUTINE radiation_init
998
999
1000!------------------------------------------------------------------------------!
1001! Description:
1002! ------------
1003!> A simple clear sky radiation model
1004!------------------------------------------------------------------------------!
1005    SUBROUTINE radiation_clearsky
1006
1007       USE indices,                                                            &
1008           ONLY:  nbgp
1009
1010       IMPLICIT NONE
1011
1012       INTEGER(iwp) :: i, j, k   !< loop indices
1013       REAL(wp)     :: exn,   &  !< Exner functions at surface
1014                       exn1,  &  !< Exner functions at first grid level
1015                       pt1       !< potential temperature at first grid level
1016
1017!
1018!--    Calculate current zenith angle
1019       CALL calc_zenith
1020
1021!
1022!--    Calculate sky transmissivity
1023       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1024
1025!
1026!--    Calculate value of the Exner function
1027       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1028!
1029!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
1030!--    point
1031       DO i = nxlg, nxrg
1032          DO j = nysg, nyng
1033             k = nzb_s_inner(j,i)
1034
1035             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1036
1037             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
1038             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1039             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1040
1041             IF ( cloud_physics )  THEN
1042                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1043                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1044             ELSE
1045                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1046             ENDIF
1047
1048             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
1049                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1050
1051          ENDDO
1052       ENDDO
1053
1054    END SUBROUTINE radiation_clearsky
1055
1056
1057!------------------------------------------------------------------------------!
1058! Description:
1059! ------------
1060!> Header output for radiation model
1061!------------------------------------------------------------------------------!
1062    SUBROUTINE radiation_header ( io )
1063
1064
1065       IMPLICIT NONE
1066 
1067       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1068   
1069
1070       
1071!
1072!--    Write radiation model header
1073       WRITE( io, 3 )
1074
1075       IF ( radiation_scheme == "constant" )  THEN
1076          WRITE( io, 4 ) net_radiation
1077       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1078          WRITE( io, 5 )
1079       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1080          WRITE( io, 6 )
1081          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1082          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1083       ENDIF
1084
1085       IF ( albedo_type == 0 )  THEN
1086          WRITE( io, 7 ) albedo
1087       ELSE
1088          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1089       ENDIF
1090       IF ( constant_albedo )  THEN
1091          WRITE( io, 9 )
1092       ENDIF
1093       
1094       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1095          WRITE ( io, 1 )  lambda
1096          WRITE ( io, 2 )  day_init, time_utc_init
1097       ENDIF
1098
1099       WRITE( io, 12 ) dt_radiation
1100 
1101
1102 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1103 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
1104            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
1105 3 FORMAT (//' Radiation model information:'/                                  &
1106              ' ----------------------------'/)
1107 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
1108           // 'W/m**2')
1109 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
1110                   ' default)')
1111 6 FORMAT ('    --> RRTMG scheme is used')
1112 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1113 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1114 9 FORMAT (/'    --> Albedo is fixed during the run')
111510 FORMAT (/'    --> Longwave radiation is disabled')
111611 FORMAT (/'    --> Shortwave radiation is disabled.')
111712 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1118
1119
1120    END SUBROUTINE radiation_header
1121   
1122
1123!------------------------------------------------------------------------------!
1124! Description:
1125! ------------
1126!> Parin for &radiation_par for radiation model
1127!------------------------------------------------------------------------------!
1128    SUBROUTINE radiation_parin
1129
1130
1131       IMPLICIT NONE
1132
1133       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1134       
1135       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
1136                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1137                                  constant_albedo, day_init, dt_radiation,     &
1138                                  lambda, lw_radiation, net_radiation,         &
1139                                  radiation_scheme, skip_time_do_radiation,    &
1140                                  sw_radiation, time_utc_init,                 &
1141                                  unscheduled_radiation_calls
1142       
1143       line = ' '
1144       
1145!
1146!--    Try to find radiation model package
1147       REWIND ( 11 )
1148       line = ' '
1149       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1150          READ ( 11, '(A)', END=10 )  line
1151       ENDDO
1152       BACKSPACE ( 11 )
1153
1154!
1155!--    Read user-defined namelist
1156       READ ( 11, radiation_par )
1157
1158!
1159!--    Set flag that indicates that the radiation model is switched on
1160       radiation = .TRUE.
1161
1162 10    CONTINUE
1163       
1164
1165    END SUBROUTINE radiation_parin
1166
1167
1168!------------------------------------------------------------------------------!
1169! Description:
1170! ------------
1171!> Implementation of the RRTMG radiation_scheme
1172!------------------------------------------------------------------------------!
1173    SUBROUTINE radiation_rrtmg
1174
1175       USE indices,                                                            &
1176           ONLY:  nbgp
1177
1178       USE particle_attributes,                                                &
1179           ONLY:  grid_particles, number_of_particles, particles,              &
1180                  particle_advection_start, prt_count
1181
1182       IMPLICIT NONE
1183
1184#if defined ( __rrtmg )
1185
1186       INTEGER(iwp) :: i, j, k, n !< loop indices
1187
1188       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
1189                        s_r3      !< weighted sum over all droplets with r^3
1190
1191!
1192!--    Calculate current (cosine of) zenith angle and whether the sun is up
1193       CALL calc_zenith     
1194!
1195!--    Calculate surface albedo
1196       IF ( .NOT. constant_albedo )  THEN
1197          CALL calc_albedo
1198       ENDIF
1199
1200!
1201!--    Prepare input data for RRTMG
1202
1203!
1204!--    In case of large scale forcing with surface data, calculate new pressure
1205!--    profile. nzt_rad might be modified by these calls and all required arrays
1206!--    will then be re-allocated
1207       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1208          CALL read_sounding_data
1209          CALL read_trace_gas_data
1210       ENDIF
1211!
1212!--    Loop over all grid points
1213       DO i = nxl, nxr
1214          DO j = nys, nyn
1215
1216!
1217!--          Prepare profiles of temperature and H2O volume mixing ratio
1218             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
1219                                                  / 1000.0_wp )**0.286_wp
1220
1221             DO k = nzb+1, nzt+1
1222                rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
1223                                 )**0.286_wp + l_d_cp * ql(k,j,i)
1224                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1225
1226             ENDDO
1227
1228!
1229!--          Avoid temperature/humidity jumps at the top of the LES domain by
1230!--          linear interpolation from nzt+2 to nzt+7
1231             DO k = nzt+2, nzt+7
1232                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
1233                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
1234                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
1235                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1236
1237                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
1238                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1239                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
1240                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1241
1242             ENDDO
1243
1244!--          Linear interpolate to zw grid
1245             DO k = nzb+2, nzt+8
1246                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
1247                                   rrtm_tlay(0,k-1))                           &
1248                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
1249                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1250             ENDDO
1251
1252
1253!
1254!--          Calculate liquid water path and cloud fraction for each column.
1255!--          Note that LWP is required in g/m² instead of kg/kg m.
1256             rrtm_cldfr  = 0.0_wp
1257             rrtm_reliq  = 0.0_wp
1258             rrtm_cliqwp = 0.0_wp
1259             rrtm_icld   = 0
1260
1261             DO k = nzb+1, nzt+1
1262                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                    &
1263                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))        &
1264                                    * 100.0_wp / g
1265
1266                IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
1267                   rrtm_cldfr(0,k) = 1.0_wp
1268                   IF ( rrtm_icld == 0 )  rrtm_icld = 1
1269
1270!
1271!--                Calculate cloud droplet effective radius
1272                   IF ( cloud_physics )  THEN
1273                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
1274                                        * rho_surface                          &
1275                                        / ( 4.0_wp * pi * nc_const * rho_l )   &
1276                                        )**0.33333333333333_wp                 &
1277                                        * EXP( LOG( sigma_gc )**2 )
1278
1279                   ELSEIF ( cloud_droplets )  THEN
1280                      number_of_particles = prt_count(k,j,i)
1281
1282                      IF (number_of_particles <= 0)  CYCLE
1283                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1284                      s_r2 = 0.0_wp
1285                      s_r3 = 0.0_wp
1286
1287                      DO  n = 1, number_of_particles
1288                         IF ( particles(n)%particle_mask )  THEN
1289                            s_r2 = s_r2 + particles(n)%radius**2 * &
1290                                   particles(n)%weight_factor
1291                            s_r3 = s_r3 + particles(n)%radius**3 * &
1292                                   particles(n)%weight_factor
1293                         ENDIF
1294                      ENDDO
1295
1296                      IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
1297
1298                   ENDIF
1299
1300!
1301!--                Limit effective radius
1302                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
1303                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1304                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1305                  ENDIF
1306                ENDIF
1307             ENDDO
1308
1309!
1310!--          Set surface temperature
1311             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1312
1313             IF ( lw_radiation )  THEN
1314               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
1315               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
1316               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
1317               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
1318               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
1319               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
1320               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
1321               rrtm_reliq      , rrtm_lw_tauaer,                               &
1322               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
1323               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
1324               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
1325
1326!
1327!--             Save fluxes
1328                DO k = nzb, nzt+1
1329                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
1330                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1331                ENDDO
1332
1333!
1334!--             Save heating rates (convert from K/d to K/h)
1335                DO k = nzb+1, nzt+1
1336                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
1337                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
1338                ENDDO
1339
1340!
1341!--             Save change in LW heating rate
1342                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1343
1344             ENDIF
1345
1346             IF ( sw_radiation .AND. sun_up )  THEN
1347                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
1348               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
1349               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
1350               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1351               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1352               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
1353               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
1354               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
1355               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
1356               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
1357               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
1358               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
1359 
1360!
1361!--             Save fluxes
1362                DO k = nzb, nzt+1
1363                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
1364                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1365                ENDDO
1366
1367!
1368!--             Save heating rates (convert from K/d to K/s)
1369                DO k = nzb+1, nzt+1
1370                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
1371                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
1372                ENDDO
1373
1374             ENDIF
1375
1376!
1377!--          Calculate surface net radiation
1378             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
1379                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1380
1381          ENDDO
1382       ENDDO
1383
1384       CALL exchange_horiz( rad_lw_in,  nbgp )
1385       CALL exchange_horiz( rad_lw_out, nbgp )
1386       CALL exchange_horiz( rad_lw_hr,    nbgp )
1387       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1388
1389       CALL exchange_horiz( rad_sw_in,  nbgp )
1390       CALL exchange_horiz( rad_sw_out, nbgp ) 
1391       CALL exchange_horiz( rad_sw_hr,    nbgp )
1392       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1393
1394       CALL exchange_horiz_2d( rad_net, nbgp )
1395       CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
1396#endif
1397
1398    END SUBROUTINE radiation_rrtmg
1399
1400
1401!------------------------------------------------------------------------------!
1402! Description:
1403! ------------
1404!> Calculate the cosine of the zenith angle (variable is called zenith)
1405!------------------------------------------------------------------------------!
1406    SUBROUTINE calc_zenith
1407
1408       IMPLICIT NONE
1409
1410       REAL(wp) ::  declination,  & !< solar declination angle
1411                    hour_angle      !< solar hour angle
1412!
1413!--    Calculate current day and time based on the initial values and simulation
1414!--    time
1415       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
1416                               / 86400.0_wp ), KIND=iwp)
1417       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1418
1419
1420!
1421!--    Calculate solar declination and hour angle   
1422       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1423       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1424
1425!
1426!--    Calculate zenith angle
1427       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
1428                                            * COS(hour_angle)
1429       zenith(0) = MAX(0.0_wp,zenith(0))
1430
1431!
1432!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
1433       IF ( zenith(0) > 0.0_wp )  THEN
1434          sun_up = .TRUE.
1435       ELSE
1436          sun_up = .FALSE.
1437       END IF
1438
1439    END SUBROUTINE calc_zenith
1440
1441#if defined ( __rrtmg ) && defined ( __netcdf )
1442!------------------------------------------------------------------------------!
1443! Description:
1444! ------------
1445!> Calculates surface albedo components based on Briegleb (1992) and
1446!> Briegleb et al. (1986)
1447!------------------------------------------------------------------------------!
1448    SUBROUTINE calc_albedo
1449
1450        IMPLICIT NONE
1451
1452        IF ( sun_up )  THEN
1453!
1454!--        Ocean
1455           IF ( albedo_type == 1 )  THEN
1456              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
1457                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
1458                                            * ( zenith(0) - 0.5_wp )           &
1459                                            * ( zenith(0) - 1.0_wp )
1460              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1461!
1462!--        Snow
1463           ELSEIF ( albedo_type == 16 )  THEN
1464              IF ( zenith(0) < 0.5_wp )  THEN
1465                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
1466                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1467                                     * zenith(0))) - 1.0_wp
1468                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
1469                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1470                                     * zenith(0))) - 1.0_wp
1471
1472                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1473                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1474              ELSE
1475                 rrtm_aldir(0,:,:) = aldif
1476                 rrtm_asdir(0,:,:) = asdif
1477              ENDIF
1478!
1479!--        Sea ice
1480           ELSEIF ( albedo_type == 15 )  THEN
1481                 rrtm_aldir(0,:,:) = aldif
1482                 rrtm_asdir(0,:,:) = asdif
1483
1484!
1485!--        Asphalt
1486           ELSEIF ( albedo_type == 17 )  THEN
1487                 rrtm_aldir(0,:,:) = aldif
1488                 rrtm_asdir(0,:,:) = asdif
1489!
1490!--        Land surfaces
1491           ELSE
1492              SELECT CASE ( albedo_type )
1493
1494!
1495!--              Surface types with strong zenith dependence
1496                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
1497                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
1498                                        (1.0_wp + 0.8_wp * zenith(0))
1499                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
1500                                        (1.0_wp + 0.8_wp * zenith(0))
1501!
1502!--              Surface types with weak zenith dependence
1503                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
1504                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
1505                                        (1.0_wp + 0.2_wp * zenith(0))
1506                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
1507                                        (1.0_wp + 0.2_wp * zenith(0))
1508
1509                 CASE DEFAULT
1510
1511              END SELECT
1512           ENDIF
1513!
1514!--        Diffusive albedo is taken from Table 2
1515           rrtm_aldif(0,:,:) = aldif
1516           rrtm_asdif(0,:,:) = asdif
1517
1518        ELSE
1519
1520           rrtm_aldir(0,:,:) = 0.0_wp
1521           rrtm_asdir(0,:,:) = 0.0_wp
1522           rrtm_aldif(0,:,:) = 0.0_wp
1523           rrtm_asdif(0,:,:) = 0.0_wp
1524        ENDIF
1525    END SUBROUTINE calc_albedo
1526
1527!------------------------------------------------------------------------------!
1528! Description:
1529! ------------
1530!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1531!------------------------------------------------------------------------------!
1532    SUBROUTINE read_sounding_data
1533
1534       IMPLICIT NONE
1535
1536       INTEGER(iwp) :: id,           & !< NetCDF id of input file
1537                       id_dim_zrad,  & !< pressure level id in the NetCDF file
1538                       id_var,       & !< NetCDF variable id
1539                       k,            & !< loop index
1540                       nz_snd,       & !< number of vertical levels in the sounding data
1541                       nz_snd_start, & !< start vertical index for sounding data to be used
1542                       nz_snd_end      !< end vertical index for souding data to be used
1543
1544       REAL(wp) :: t_surface           !< actual surface temperature
1545
1546       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1547                                               t_snd_tmp      !< temporary temperature profile (sounding)
1548
1549!
1550!--    In case of updates, deallocate arrays first (sufficient to check one
1551!--    array as the others are automatically allocated). This is required
1552!--    because nzt_rad might change during the update
1553       IF ( ALLOCATED ( hyp_snd ) )  THEN
1554          DEALLOCATE( hyp_snd )
1555          DEALLOCATE( t_snd )
1556          DEALLOCATE( q_snd  )
1557          DEALLOCATE ( rrtm_play )
1558          DEALLOCATE ( rrtm_plev )
1559          DEALLOCATE ( rrtm_tlay )
1560          DEALLOCATE ( rrtm_tlev )
1561
1562          DEALLOCATE ( rrtm_h2ovmr )
1563          DEALLOCATE ( rrtm_cicewp )
1564          DEALLOCATE ( rrtm_cldfr )
1565          DEALLOCATE ( rrtm_cliqwp )
1566          DEALLOCATE ( rrtm_reice )
1567          DEALLOCATE ( rrtm_reliq )
1568          DEALLOCATE ( rrtm_lw_taucld )
1569          DEALLOCATE ( rrtm_lw_tauaer )
1570
1571          DEALLOCATE ( rrtm_lwdflx  )
1572          DEALLOCATE ( rrtm_lwdflxc )
1573          DEALLOCATE ( rrtm_lwuflx  )
1574          DEALLOCATE ( rrtm_lwuflxc )
1575          DEALLOCATE ( rrtm_lwuflx_dt )
1576          DEALLOCATE ( rrtm_lwuflxc_dt )
1577          DEALLOCATE ( rrtm_lwhr  )
1578          DEALLOCATE ( rrtm_lwhrc )
1579
1580          DEALLOCATE ( rrtm_sw_taucld )
1581          DEALLOCATE ( rrtm_sw_ssacld )
1582          DEALLOCATE ( rrtm_sw_asmcld )
1583          DEALLOCATE ( rrtm_sw_fsfcld )
1584          DEALLOCATE ( rrtm_sw_tauaer )
1585          DEALLOCATE ( rrtm_sw_ssaaer )
1586          DEALLOCATE ( rrtm_sw_asmaer ) 
1587          DEALLOCATE ( rrtm_sw_ecaer )   
1588 
1589          DEALLOCATE ( rrtm_swdflx  )
1590          DEALLOCATE ( rrtm_swdflxc )
1591          DEALLOCATE ( rrtm_swuflx  )
1592          DEALLOCATE ( rrtm_swuflxc )
1593          DEALLOCATE ( rrtm_swhr  )
1594          DEALLOCATE ( rrtm_swhrc )
1595
1596       ENDIF
1597
1598!
1599!--    Open file for reading
1600       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1601       CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1602
1603!
1604!--    Inquire dimension of z axis and save in nz_snd
1605       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1606       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1607       CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1608
1609!
1610! !--    Allocate temporary array for storing pressure data
1611       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1612       hyp_snd_tmp = 0.0_wp
1613
1614
1615!--    Read pressure from file
1616       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1617       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
1618                               count = (/nz_snd/) )
1619       CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1620
1621!
1622!--    Allocate temporary array for storing temperature data
1623       ALLOCATE( t_snd_tmp(1:nz_snd) )
1624       t_snd_tmp = 0.0_wp
1625
1626!
1627!--    Read temperature from file
1628       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1629       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
1630                               count = (/nz_snd/) )
1631       CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1632
1633!
1634!--    Calculate start of sounding data
1635       nz_snd_start = nz_snd + 1
1636       nz_snd_end   = nz_snd + 1
1637
1638!
1639!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
1640!--    in Pa, hyp_snd in hPa).
1641       DO  k = 1, nz_snd
1642          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
1643             nz_snd_start = k
1644             EXIT
1645          END IF
1646       END DO
1647
1648       IF ( nz_snd_start <= nz_snd )  THEN
1649          nz_snd_end = nz_snd
1650       END IF
1651
1652
1653!
1654!--    Calculate of total grid points for RRTMG calculations
1655       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1656
1657!
1658!--    Save data above LES domain in hyp_snd, t_snd and q_snd
1659!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1660       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1661       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
1662       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
1663       hyp_snd = 0.0_wp
1664       t_snd = 0.0_wp
1665       q_snd = 0.0_wp
1666
1667       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1668       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1669
1670       nc_stat = NF90_CLOSE( id )
1671
1672!
1673!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
1674!--    top of the LES domain. This routine does not consider horizontal or
1675!--    vertical variability of pressure and temperature
1676       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
1677       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
1678
1679       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1680       DO k = nzb+1, nzt+1
1681          rrtm_play(0,k) = hyp(k) * 0.01_wp
1682          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1683                         t_surface )**(1.0_wp/0.286_wp)
1684       ENDDO
1685
1686       DO k = nzt+2, nzt_rad
1687          rrtm_play(0,k) = hyp_snd(k)
1688          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1689       ENDDO
1690       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
1691                                   1.5 * hyp_snd(nzt_rad)                      &
1692                                 - 0.5 * hyp_snd(nzt_rad-1) )
1693       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
1694                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
1695
1696       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1697
1698!
1699!--    Calculate temperature/humidity levels at top of the LES domain.
1700!--    Currently, the temperature is taken from sounding data (might lead to a
1701!--    temperature jump at interface. To do: Humidity is currently not
1702!--    calculated above the LES domain.
1703       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
1704       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
1705       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1706
1707       DO k = nzt+8, nzt_rad
1708          rrtm_tlay(0,k)   = t_snd(k)
1709          rrtm_h2ovmr(0,k) = q_snd(k)
1710       ENDDO
1711       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
1712                                - rrtm_tlay(0,nzt_rad-1)
1713       DO k = nzt+9, nzt_rad+1
1714          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
1715                             - rrtm_tlay(0,k-1))                               &
1716                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
1717                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1718       ENDDO
1719       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1720
1721       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
1722                                  - rrtm_tlev(0,nzt_rad)
1723!
1724!--    Allocate remaining RRTMG arrays
1725       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1726       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1727       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1728       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1729       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1730       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1731       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1732       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1733       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1734       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1735       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1736       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1737       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1738       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
1739       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
1740
1741!
1742!--    The ice phase is currently not considered in PALM
1743       rrtm_cicewp = 0.0_wp
1744       rrtm_reice  = 0.0_wp
1745
1746!
1747!--    Set other parameters (move to NAMELIST parameters in the future)
1748       rrtm_lw_tauaer = 0.0_wp
1749       rrtm_lw_taucld = 0.0_wp
1750       rrtm_sw_taucld = 0.0_wp
1751       rrtm_sw_ssacld = 0.0_wp
1752       rrtm_sw_asmcld = 0.0_wp
1753       rrtm_sw_fsfcld = 0.0_wp
1754       rrtm_sw_tauaer = 0.0_wp
1755       rrtm_sw_ssaaer = 0.0_wp
1756       rrtm_sw_asmaer = 0.0_wp
1757       rrtm_sw_ecaer  = 0.0_wp
1758
1759
1760       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
1761       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
1762       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
1763       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1764       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1765       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1766
1767       rrtm_swdflx  = 0.0_wp
1768       rrtm_swuflx  = 0.0_wp
1769       rrtm_swhr    = 0.0_wp 
1770       rrtm_swuflxc = 0.0_wp
1771       rrtm_swdflxc = 0.0_wp
1772       rrtm_swhrc   = 0.0_wp
1773
1774       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
1775       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
1776       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
1777       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1778       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1779       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1780
1781       rrtm_lwdflx  = 0.0_wp
1782       rrtm_lwuflx  = 0.0_wp
1783       rrtm_lwhr    = 0.0_wp 
1784       rrtm_lwuflxc = 0.0_wp
1785       rrtm_lwdflxc = 0.0_wp
1786       rrtm_lwhrc   = 0.0_wp
1787
1788       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
1789       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
1790
1791       rrtm_lwuflx_dt = 0.0_wp
1792       rrtm_lwuflxc_dt = 0.0_wp
1793
1794    END SUBROUTINE read_sounding_data
1795
1796
1797!------------------------------------------------------------------------------!
1798! Description:
1799! ------------
1800!> Read trace gas data from file
1801!------------------------------------------------------------------------------!
1802    SUBROUTINE read_trace_gas_data
1803
1804       USE rrsw_ncpar
1805
1806       IMPLICIT NONE
1807
1808       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
1809
1810       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
1811           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
1812                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
1813
1814       INTEGER(iwp) :: id,     & !< NetCDF id
1815                       k,      & !< loop index
1816                       m,      & !< loop index
1817                       n,      & !< loop index
1818                       nabs,   & !< number of absorbers
1819                       np,     & !< number of pressure levels
1820                       id_abs, & !< NetCDF id of the respective absorber
1821                       id_dim, & !< NetCDF id of asborber's dimension
1822                       id_var    !< NetCDf id ot the absorber
1823
1824       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
1825
1826
1827       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
1828                                                 rrtm_play_tmp, & !< temporary array for pressure zu-levels
1829                                                 rrtm_plev_tmp, & !< temporary array for pressure zw-levels
1830                                                 trace_path_tmp   !< temporary array for storing trace gas path data
1831
1832       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
1833                                                 trace_mls_path, & !< array for storing trace gas path data
1834                                                 trace_mls_tmp     !< temporary array for storing trace gas data
1835
1836
1837!
1838!--    In case of updates, deallocate arrays first (sufficient to check one
1839!--    array as the others are automatically allocated)
1840       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
1841          DEALLOCATE ( rrtm_o3vmr  )
1842          DEALLOCATE ( rrtm_co2vmr )
1843          DEALLOCATE ( rrtm_ch4vmr )
1844          DEALLOCATE ( rrtm_n2ovmr )
1845          DEALLOCATE ( rrtm_o2vmr  )
1846          DEALLOCATE ( rrtm_cfc11vmr )
1847          DEALLOCATE ( rrtm_cfc12vmr )
1848          DEALLOCATE ( rrtm_cfc22vmr )
1849          DEALLOCATE ( rrtm_ccl4vmr  )
1850       ENDIF
1851
1852!
1853!--    Allocate trace gas profiles
1854       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
1855       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
1856       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
1857       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
1858       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
1859       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
1860       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
1861       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
1862       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
1863
1864!
1865!--    Open file for reading
1866       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1867       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
1868!
1869!--    Inquire dimension ids and dimensions
1870       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
1871       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1872       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
1873       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1874
1875       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
1876       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1877       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
1878       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1879   
1880
1881!
1882!--    Allocate pressure, and trace gas arrays     
1883       ALLOCATE( p_mls(1:np) )
1884       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
1885       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
1886
1887
1888       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1889       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1890       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
1891       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1892
1893       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
1894       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1895       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
1896       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1897
1898
1899!
1900!--    Write absorber amounts (mls) to trace_mls
1901       DO n = 1, num_trace_gases
1902          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
1903
1904          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
1905
1906!
1907!--       Replace missing values by zero
1908          WHERE ( trace_mls(n,:) > 2.0_wp ) 
1909             trace_mls(n,:) = 0.0_wp
1910          END WHERE
1911       END DO
1912
1913       DEALLOCATE ( trace_mls_tmp )
1914
1915       nc_stat = NF90_CLOSE( id )
1916       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
1917
1918!
1919!--    Add extra pressure level for calculations of the trace gas paths
1920       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
1921       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
1922
1923       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
1924       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
1925       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
1926       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
1927                                         * rrtm_plev(0,nzt_rad+1) )
1928 
1929!
1930!--    Calculate trace gas path (zero at surface) with interpolation to the
1931!--    sounding levels
1932       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
1933
1934       trace_mls_path(nzb+1,:) = 0.0_wp
1935       
1936       DO k = nzb+2, nzt_rad+2
1937          DO m = 1, num_trace_gases
1938             trace_mls_path(k,m) = trace_mls_path(k-1,m)
1939
1940!
1941!--          When the pressure level is higher than the trace gas pressure
1942!--          level, assume that
1943             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
1944               
1945                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
1946                                      * ( rrtm_plev_tmp(k-1)                   &
1947                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
1948                                        ) / g
1949             ENDIF
1950
1951!
1952!--          Integrate for each sounding level from the contributing p_mls
1953!--          levels
1954             DO n = 2, np
1955!
1956!--             Limit p_mls so that it is within the model level
1957                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
1958                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
1959                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
1960                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
1961
1962                IF ( p_mls_l > p_mls_u )  THEN
1963
1964!
1965!--                Calculate weights for interpolation
1966                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
1967                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
1968                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
1969
1970!
1971!--                Add level to trace gas path
1972                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
1973                                         +  ( p_wgt_u * trace_mls(m,n)         &
1974                                            + p_wgt_l * trace_mls(m,n-1) )     &
1975                                         * (p_mls_l - p_mls_u) / g
1976                ENDIF
1977             ENDDO
1978
1979             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
1980                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
1981                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
1982                                          - rrtm_plev_tmp(k)                   &
1983                                        ) / g
1984             ENDIF 
1985          ENDDO
1986       ENDDO
1987
1988
1989!
1990!--    Prepare trace gas path profiles
1991       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
1992
1993       DO m = 1, num_trace_gases
1994
1995          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
1996                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
1997                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
1998                                       - rrtm_plev_tmp(2:nzt_rad+2) )
1999
2000!
2001!--       Save trace gas paths to the respective arrays
2002          SELECT CASE ( TRIM( trace_names(m) ) )
2003
2004             CASE ( 'O3' )
2005
2006                rrtm_o3vmr(0,:) = trace_path_tmp(:)
2007
2008             CASE ( 'CO2' )
2009
2010                rrtm_co2vmr(0,:) = trace_path_tmp(:)
2011
2012             CASE ( 'CH4' )
2013
2014                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2015
2016             CASE ( 'N2O' )
2017
2018                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2019
2020             CASE ( 'O2' )
2021
2022                rrtm_o2vmr(0,:) = trace_path_tmp(:)
2023
2024             CASE ( 'CFC11' )
2025
2026                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2027
2028             CASE ( 'CFC12' )
2029
2030                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2031
2032             CASE ( 'CFC22' )
2033
2034                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2035
2036             CASE ( 'CCL4' )
2037
2038                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2039
2040             CASE DEFAULT
2041
2042          END SELECT
2043
2044       ENDDO
2045
2046       DEALLOCATE ( trace_path_tmp )
2047       DEALLOCATE ( trace_mls_path )
2048       DEALLOCATE ( rrtm_play_tmp )
2049       DEALLOCATE ( rrtm_plev_tmp )
2050       DEALLOCATE ( trace_mls )
2051       DEALLOCATE ( p_mls )
2052
2053    END SUBROUTINE read_trace_gas_data
2054
2055
2056    SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2057
2058       USE control_parameters,                                                 &
2059           ONLY:  message_string
2060
2061       USE NETCDF
2062
2063       USE pegrid
2064
2065       IMPLICIT NONE
2066
2067       CHARACTER(LEN=6) ::  message_identifier
2068       CHARACTER(LEN=*) ::  routine_name
2069
2070       INTEGER(iwp) ::  errno
2071
2072       IF ( nc_stat /= NF90_NOERR )  THEN
2073
2074          WRITE( message_identifier, '(''NC'',I4.4)' )  errno
2075          message_string = TRIM( NF90_STRERROR( nc_stat ) )
2076
2077          CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2078
2079       ENDIF
2080
2081    END SUBROUTINE netcdf_handle_error_rad
2082#endif
2083
2084
2085!------------------------------------------------------------------------------!
2086! Description:
2087! ------------
2088!> Calculate temperature tendency due to radiative cooling/heating.
2089!> Cache-optimized version.
2090!------------------------------------------------------------------------------!
2091    SUBROUTINE radiation_tendency_ij ( i, j, tend )
2092
2093       USE cloud_parameters,                                                   &
2094           ONLY:  pt_d_t
2095
2096       IMPLICIT NONE
2097
2098       INTEGER(iwp) :: i, j, k !< loop indices
2099
2100       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2101
2102#if defined ( __rrtmg )
2103!
2104!--    Calculate tendency based on heating rate
2105       DO k = nzb+1, nzt+1
2106          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
2107                                      * pt_d_t(k) * d_seconds_hour
2108       ENDDO
2109
2110#endif
2111
2112    END SUBROUTINE radiation_tendency_ij
2113
2114
2115!------------------------------------------------------------------------------!
2116! Description:
2117! ------------
2118!> Calculate temperature tendency due to radiative cooling/heating.
2119!> Vector-optimized version
2120!------------------------------------------------------------------------------!
2121    SUBROUTINE radiation_tendency ( tend )
2122
2123       USE cloud_parameters,                                                   &
2124           ONLY:  pt_d_t
2125
2126       USE indices,                                                            &
2127           ONLY:  nxl, nxr, nyn, nys
2128
2129       IMPLICIT NONE
2130
2131       INTEGER(iwp) :: i, j, k !< loop indices
2132
2133       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2134
2135#if defined ( __rrtmg )
2136!
2137!--    Calculate tendency based on heating rate
2138       DO  i = nxl, nxr
2139          DO  j = nys, nyn
2140             DO k = nzb+1, nzt+1
2141                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
2142                                            +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
2143                                            * d_seconds_hour
2144             ENDDO
2145         ENDDO
2146       ENDDO
2147#endif
2148
2149    END SUBROUTINE radiation_tendency
2150
2151 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.