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

Last change on this file since 1856 was 1856, checked in by maronga, 5 years ago

bugfix in land surface / radiation model

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