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

Last change on this file since 1952 was 1857, checked in by maronga, 8 years ago

last commit documented

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