source: palm/trunk/SOURCE/radiation_model.f90 @ 1702

Last change on this file since 1702 was 1702, checked in by maronga, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 68.3 KB
Line 
1!> @file radiation_model.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-2015 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: radiation_model.f90 1702 2015-11-02 07:44:12Z maronga $
26!
27! 1701 2015-11-02 07:43:04Z maronga
28! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
29!
30! 1691 2015-10-26 16:17:44Z maronga
31! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
32! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
33! Added output of radiative heating rates.
34!
35! 1682 2015-10-07 23:56:08Z knoop
36! Code annotations made doxygen readable
37!
38! 1606 2015-06-29 10:43:37Z maronga
39! Added preprocessor directive __netcdf to allow for compiling without netCDF.
40! Note, however, that RRTMG cannot be used without netCDF.
41!
42! 1590 2015-05-08 13:56:27Z maronga
43! Bugfix: definition of character strings requires same length for all elements
44!
45! 1587 2015-05-04 14:19:01Z maronga
46! Added albedo class for snow
47!
48! 1585 2015-04-30 07:05:52Z maronga
49! Added support for RRTMG
50!
51! 1571 2015-03-12 16:12:49Z maronga
52! Added missing KIND attribute. Removed upper-case variable names
53!
54! 1551 2015-03-03 14:18:16Z maronga
55! Added support for data output. Various variables have been renamed. Added
56! interface for different radiation schemes (currently: clear-sky, constant, and
57! RRTM (not yet implemented).
58!
59! 1496 2014-12-02 17:25:50Z maronga
60! Initial revision
61!
62!
63! Description:
64! ------------
65!> Radiation models and interfaces
66!> @todo move variable definitions used in init_radiation only to the subroutine
67!>       as they are no longer required after initialization.
68!> @todo Output of full column vertical profiles used in RRTMG
69!> @todo Output of other rrtm arrays (such as volume mixing ratios)
70!> @todo Adapt for use with topography
71!>
72!> @note Many variables have a leading dummy dimension (0:0) in order to
73!>       match the assume-size shape expected by the RRTMG model.
74!------------------------------------------------------------------------------!
75 MODULE radiation_model_mod
76 
77
78    USE arrays_3d,                                                             &
79        ONLY:  dzw, hyp, pt, q, ql, zw
80
81    USE cloud_parameters,                                                      &
82        ONLY:  cp, l_d_cp, nc_const, rho_l, sigma_gc 
83
84    USE constants,                                                             &
85        ONLY:  pi
86
87    USE control_parameters,                                                    &
88        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
89               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
90               surface_pressure, time_since_reference_point
91
92    USE indices,                                                               &
93        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
94
95    USE kinds
96
97#if defined ( __netcdf )
98    USE netcdf
99#endif
100
101    USE netcdf_control,                                                        &
102        ONLY:  dots_label, dots_num, dots_unit
103
104#if defined ( __rrtmg )
105    USE parrrsw,                                                               &
106        ONLY:  naerec, nbndsw
107
108    USE parrrtm,                                                               &
109        ONLY:  nbndlw
110
111    USE rrtmg_lw_init,                                                         &
112        ONLY:  rrtmg_lw_ini
113
114    USE rrtmg_sw_init,                                                         &
115        ONLY:  rrtmg_sw_ini
116
117    USE rrtmg_lw_rad,                                                          &
118        ONLY:  rrtmg_lw
119
120    USE rrtmg_sw_rad,                                                          &
121        ONLY:  rrtmg_sw
122#endif
123
124
125
126    IMPLICIT NONE
127
128    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
129
130!
131!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
132    CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/      &
133                                   'user defined                         ', & !  0
134                                   'ocean                                ', & !  1
135                                   'mixed farming, tall grassland        ', & !  2
136                                   'tall/medium grassland                ', & !  3
137                                   'evergreen shrubland                  ', & !  4
138                                   'short grassland/meadow/shrubland     ', & !  5
139                                   'evergreen needleleaf forest          ', & !  6
140                                   'mixed deciduous evergreen forest     ', & !  7
141                                   'deciduous forest                     ', & !  8
142                                   'tropical evergreen broadleaved forest', & !  9
143                                   'medium/tall grassland/woodland       ', & ! 10
144                                   'desert, sandy                        ', & ! 11
145                                   'desert, rocky                        ', & ! 12
146                                   'tundra                               ', & ! 13
147                                   'land ice                             ', & ! 14
148                                   'sea ice                              ', & ! 15
149                                   'snow                                 '  & ! 16
150                                                         /)
151
152    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
153                    day,                 & !< current day of the year
154                    day_init     = 172,  & !< day of the year at model start (21/06)
155                    dots_rad     = 0       !< starting index for timeseries output
156
157
158
159
160
161
162    LOGICAL ::  constant_albedo = .FALSE.,       & !< flag parameter indicating whether the albedo may change depending on zenith
163                force_radiation_call = .FALSE.,  & !< flag parameter for unscheduled radiation calls
164                lw_radiation = .TRUE.,           & !< flag parameter indicing whether longwave radiation shall be calculated
165                radiation = .FALSE.,             & !< flag parameter indicating whether the radiation model is used
166                sun_up    = .TRUE.,              & !< flag parameter indicating whether the sun is up or down
167                sw_radiation = .TRUE.              !< flag parameter indicing whether shortwave radiation shall be calculated
168
169
170    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
171                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
172                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
173                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
174
175    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
176                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
177                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
178                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
179                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
180                decl_1,                          & !< declination coef. 1
181                decl_2,                          & !< declination coef. 2
182                decl_3,                          & !< declination coef. 3
183                dt_radiation = 0.0_wp,           & !< radiation model timestep
184                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
185                lambda = 0.0_wp,                 & !< longitude in degrees
186                lon = 0.0_wp,                    & !< longitude in radians
187                lat = 0.0_wp,                    & !< latitude in radians
188                net_radiation = 0.0_wp,          & !< net radiation at surface
189                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
190                sky_trans,                       & !< sky transmissivity
191                time_radiation = 0.0_wp,         & !< time since last call of radiation code
192                time_utc,                        & !< current time in UTC
193                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
194
195    REAL(wp), DIMENSION(0:0) ::  zenith        !< solar zenith angle
196
197    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
198                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
199                rad_net,                     & !< net radiation at the surface
200                rad_net_av                     !< average of rad_net
201
202!
203!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
204!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
205    REAL(wp), DIMENSION(0:2,1:16), PARAMETER :: albedo_pars = RESHAPE( (/& 
206                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
207                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
208                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
209                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
210                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
211                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
212                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
213                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
214                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
215                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
216                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
217                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
218                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
219                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
220                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
221                                   0.95_wp, 0.70_wp, 0.82_wp             & ! 16
222                                 /), (/ 3, 16 /) )
223
224    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
225                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
226                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
227                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
228                        rad_lw_hr_av,                  & !< average of rad_sw_hr
229                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
230                        rad_lw_in_av,                  & !< average of rad_lw_in
231                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
232                        rad_lw_out_av,                 & !< average of rad_lw_out
233                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
234                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
235                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
236                        rad_sw_hr_av,                  & !< average of rad_sw_hr
237                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
238                        rad_sw_in_av,                  & !< average of rad_sw_in
239                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
240                        rad_sw_out_av                    !< average of rad_sw_out
241
242
243!
244!-- Variables and parameters used in RRTMG only
245#if defined ( __rrtmg )
246    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
247
248
249!
250!-- Flag parameters for RRTMGS (should not be changed)
251    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
252                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
253                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
254                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
255                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
256                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
257
258!
259!-- The following variables should be only changed with care, as this will
260!-- require further setting of some variables, which is currently not
261!-- implemented (aerosols, ice phase).
262    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
263                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
264                    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)
265                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
266
267    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
268
269    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
270
271    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
272                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
273                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
274                                           t_snd          !< actual temperature from sounding data (hPa)
275
276    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
277                                             aldir,          & !< longwave direct albedo solar angle of 60°
278                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
279                                             asdir,          & !< shortwave direct albedo solar angle of 60°
280                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
281                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
282                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
283                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
284                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
285                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
286                                             rrtm_cldfr,     & !< cloud fraction (0,1)
287                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
288                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
289                                             rrtm_emis,      & !< surface emissivity (0-1)   
290                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
291                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
292                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
293                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
294                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
295                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
296                                             rrtm_reice,     & !< cloud ice effective radius (microns)
297                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
298                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
299                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
300                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
301                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
302                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
303                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
304                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
305                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
306                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
307                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
308                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
309                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
310                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
311                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
312                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
313                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
314
315!
316!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
317    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
318                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
319                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
320                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
321                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
322                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
323                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
324                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
325                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
326                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
327                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
328                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
329                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
330                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
331                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
332                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
333                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
334                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
335
336#endif
337
338    INTERFACE init_radiation
339       MODULE PROCEDURE init_radiation
340    END INTERFACE init_radiation
341
342    INTERFACE radiation_clearsky
343       MODULE PROCEDURE radiation_clearsky
344    END INTERFACE radiation_clearsky
345
346    INTERFACE radiation_rrtmg
347       MODULE PROCEDURE radiation_rrtmg
348    END INTERFACE radiation_rrtmg
349
350    INTERFACE radiation_tendency
351       MODULE PROCEDURE radiation_tendency
352       MODULE PROCEDURE radiation_tendency_ij
353    END INTERFACE radiation_tendency
354
355    SAVE
356
357    PRIVATE
358
359    PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
360           albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad,  &
361           dt_radiation, emissivity, force_radiation_call, init_radiation,     &
362           lambda, lw_radiation, net_radiation, rad_net, rad_net_av, radiation,&
363           radiation_clearsky, radiation_rrtmg, radiation_scheme,              &
364           radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out,            &
365           rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,            &
366           rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,   &
367           rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,   &
368           skip_time_do_radiation, sw_radiation, time_radiation, time_utc_init
369
370
371#if defined ( __rrtmg )
372    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv,          &
373           rrtm_lwuflx_dt
374#endif
375
376 CONTAINS
377
378!------------------------------------------------------------------------------!
379! Description:
380! ------------
381!> Initialization of the radiation model
382!------------------------------------------------------------------------------!
383    SUBROUTINE init_radiation
384   
385       IMPLICIT NONE
386
387!
388!--    Allocate array for storing the surface net radiation
389       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
390          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
391          rad_net = 0.0_wp
392       ENDIF
393
394!
395!--    Fix net radiation in case of radiation_scheme = 'constant'
396       IF ( radiation_scheme == 'constant' )  THEN
397          rad_net = net_radiation
398          radiation = .FALSE.
399!
400!--    Calculate orbital constants
401       ELSE
402          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
403          decl_2 = 2.0_wp * pi / 365.0_wp
404          decl_3 = decl_2 * 81.0_wp
405          lat    = phi * pi / 180.0_wp
406          lon    = lambda * pi / 180.0_wp
407       ENDIF
408
409
410       IF ( radiation_scheme == 'clear-sky' )  THEN
411
412          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
413
414          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
415             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
416          ENDIF
417          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
418             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
419          ENDIF
420
421          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
422             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
423          ENDIF
424          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
425             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
426          ENDIF
427
428          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
429             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
430          ENDIF
431          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
432             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
433          ENDIF
434
435          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
436             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
437          ENDIF
438          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
439             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
440          ENDIF
441
442          rad_sw_in  = 0.0_wp
443          rad_sw_out = 0.0_wp
444          rad_lw_in  = 0.0_wp
445          rad_lw_out = 0.0_wp
446
447!
448!--       Overwrite albedo if manually set in parameter file
449          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
450             albedo = albedo_pars(2,albedo_type)
451          ENDIF
452   
453          alpha = albedo
454 
455!
456!--    Initialization actions for RRTMG
457       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
458#if defined ( __rrtmg )
459!
460!--       Allocate albedos
461          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
462          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
463          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
464          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
465          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
466          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
467          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
468          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
469
470          IF ( albedo_type /= 0 )  THEN
471             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
472                albedo_lw_dif = albedo_pars(0,albedo_type)
473                albedo_lw_dir = albedo_lw_dif
474             ENDIF
475             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
476                albedo_sw_dif = albedo_pars(1,albedo_type)
477                albedo_sw_dir = albedo_sw_dif
478             ENDIF
479          ENDIF
480
481          aldif(:,:) = albedo_lw_dif
482          aldir(:,:) = albedo_lw_dir
483          asdif(:,:) = albedo_sw_dif
484          asdir(:,:) = albedo_sw_dir
485!
486!--       Calculate initial values of current (cosine of) the zenith angle and
487!--       whether the sun is up
488          CALL calc_zenith     
489!
490!--       Calculate initial surface albedo
491          IF ( .NOT. constant_albedo )  THEN
492             CALL calc_albedo
493          ELSE
494             rrtm_aldif(0,:,:) = aldif(:,:)
495             rrtm_aldir(0,:,:) = aldir(:,:)
496             rrtm_asdif(0,:,:) = asdif(:,:) 
497             rrtm_asdir(0,:,:) = asdir(:,:)   
498          ENDIF
499
500!
501!--       Allocate surface emissivity
502          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
503          rrtm_emis = emissivity
504
505!
506!--       Allocate 3d arrays of radiative fluxes and heating rates
507          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
508             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
509             rad_sw_in = 0.0_wp
510          ENDIF
511
512          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
513             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
514          ENDIF
515
516          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
517             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
518             rad_sw_out = 0.0_wp
519          ENDIF
520
521          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
522             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
523          ENDIF
524
525          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
526             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
527             rad_sw_hr = 0.0_wp
528          ENDIF
529
530          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
531             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
532             rad_sw_hr_av = 0.0_wp
533          ENDIF
534
535          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
536             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
537             rad_sw_cs_hr = 0.0_wp
538          ENDIF
539
540          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
541             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
542             rad_sw_cs_hr_av = 0.0_wp
543          ENDIF
544
545          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
546             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
547             rad_lw_in     = 0.0_wp
548          ENDIF
549
550          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
551             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
552          ENDIF
553
554          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
555             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
556            rad_lw_out    = 0.0_wp
557          ENDIF
558
559          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
560             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
561          ENDIF
562
563          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
564             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
565             rad_lw_hr = 0.0_wp
566          ENDIF
567
568          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
569             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
570             rad_lw_hr_av = 0.0_wp
571          ENDIF
572
573          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
574             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
575             rad_lw_cs_hr = 0.0_wp
576          ENDIF
577
578          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
579             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
580             rad_lw_cs_hr_av = 0.0_wp
581          ENDIF
582
583          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
584          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
585          rad_sw_cs_in  = 0.0_wp
586          rad_sw_cs_out = 0.0_wp
587
588          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
589          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
590          rad_lw_cs_in  = 0.0_wp
591          rad_lw_cs_out = 0.0_wp
592
593!
594!--       Allocate dummy array for storing surface temperature
595          ALLOCATE ( rrtm_tsfc(1) )
596
597!
598!--       Initialize RRTMG
599          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
600          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
601
602!
603!--       Set input files for RRTMG
604          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
605          IF ( .NOT. snd_exists )  THEN
606             rrtm_input_file = "rrtmg_lw.nc"
607          ENDIF
608
609!
610!--       Read vertical layers for RRTMG from sounding data
611!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
612!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
613!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
614          CALL read_sounding_data
615
616!
617!--       Read trace gas profiles from file. This routine provides
618!--       the rrtm_ arrays (1:nzt_rad+1)
619          CALL read_trace_gas_data
620#endif
621       ENDIF
622
623!
624!--    Perform user actions if required
625       CALL user_init_radiation
626
627
628!
629!--    Add timeseries for radiation model
630       dots_rad = dots_num + 1
631       dots_num = dots_num + 5
632
633       dots_label(dots_rad) = "rad_net"
634       dots_label(dots_rad+1) = "rad_lw_in"
635       dots_label(dots_rad+2) = "rad_lw_out"
636       dots_label(dots_rad+3) = "rad_sw_in"
637       dots_label(dots_rad+4) = "rad_sw_out"
638       dots_unit(dots_rad:dots_rad+4) = "W/m2"
639
640!
641!--    Output of albedos is only required for RRTMG
642       IF ( radiation_scheme == 'rrtmg' )  THEN
643          dots_num  = dots_num + 4
644          dots_label(dots_rad+5) = "rrtm_aldif"
645          dots_label(dots_rad+6) = "rrtm_aldir"
646          dots_label(dots_rad+7) = "rrtm_asdif"
647          dots_label(dots_rad+8) = "rrtm_asdir"
648          dots_unit(dots_num+5:dots_num+8) = ""
649
650       ENDIF
651
652!
653!--    Calculate radiative fluxes at model start
654       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
655          IF ( radiation_scheme == 'clear-sky' )  THEN
656              CALL radiation_clearsky
657          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
658             CALL radiation_rrtmg
659          ENDIF
660       ENDIF
661
662       RETURN
663
664    END SUBROUTINE init_radiation
665
666
667!------------------------------------------------------------------------------!
668! Description:
669! ------------
670!> A simple clear sky radiation model
671!------------------------------------------------------------------------------!
672    SUBROUTINE radiation_clearsky
673
674       USE indices,                                                            &
675           ONLY:  nbgp
676
677       IMPLICIT NONE
678
679       INTEGER(iwp) :: i, j, k   !< loop indices
680       REAL(wp)     :: exn,   &  !< Exner functions at surface
681                       exn_1, &  !< Exner functions at first grid level
682                       pt_1      !< potential temperature at first grid level
683
684!
685!--    Calculate current zenith angle
686       CALL calc_zenith
687
688!
689!--    Calculate sky transmissivity
690       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
691
692!
693!--    Calculate value of the Exner function
694       exn = (surface_pressure / 1000.0_wp )**0.286_wp
695!
696!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
697!--    point
698       DO i = nxl, nxr
699          DO j = nys, nyn
700             k = nzb_s_inner(j,i)
701
702             exn_1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
703
704             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
705             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
706             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
707
708             IF ( cloud_physics )  THEN
709                pt_1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)
710                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt_1 * exn_1)**4
711             ELSE
712                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn_1)**4
713             ENDIF
714
715             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
716                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
717
718          ENDDO
719       ENDDO
720
721       CALL exchange_horiz_2d( rad_lw_in,  nbgp )
722       CALL exchange_horiz_2d( rad_lw_out, nbgp )
723       CALL exchange_horiz_2d( rad_sw_in,  nbgp )
724       CALL exchange_horiz_2d( rad_sw_out, nbgp )
725       CALL exchange_horiz_2d( rad_net,    nbgp )
726
727       RETURN
728
729    END SUBROUTINE radiation_clearsky
730
731
732!------------------------------------------------------------------------------!
733! Description:
734! ------------
735!> Implementation of the RRTMG radiation_scheme
736!------------------------------------------------------------------------------!
737    SUBROUTINE radiation_rrtmg
738
739       USE indices,                                                            &
740           ONLY:  nbgp
741
742       USE particle_attributes,                                                &
743           ONLY:  grid_particles, number_of_particles, particles,              &
744                  particle_advection_start, prt_count
745
746       IMPLICIT NONE
747
748#if defined ( __rrtmg )
749
750       INTEGER(iwp) :: i, j, k, n !< loop indices
751
752       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
753                        s_r3      !< weighted sum over all droplets with r^3
754
755!
756!--    Calculate current (cosine of) zenith angle and whether the sun is up
757       CALL calc_zenith     
758!
759!--    Calculate surface albedo
760       IF ( .NOT. constant_albedo )  THEN
761          CALL calc_albedo
762       ENDIF
763
764!
765!--    Prepare input data for RRTMG
766
767!
768!--    In case of large scale forcing with surface data, calculate new pressure
769!--    profile. nzt_rad might be modified by these calls and all required arrays
770!--    will then be re-allocated
771       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
772          CALL read_sounding_data
773          CALL read_trace_gas_data
774       ENDIF
775!
776!--    Loop over all grid points
777       DO i = nxl, nxr
778          DO j = nys, nyn
779
780!
781!--          Prepare profiles of temperature and H2O volume mixing ratio
782             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
783                                                  / 1000.0_wp )**0.286_wp
784
785             DO k = nzb+1, nzt+1
786                rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
787                                 )**0.286_wp + l_d_cp * ql(k,j,i)
788                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
789
790             ENDDO
791
792!
793!--          Avoid temperature/humidity jumps at the top of the LES domain by
794!--          linear interpolation from nzt+2 to nzt+7
795             DO k = nzt+2, nzt+7
796                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
797                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
798                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
799                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
800
801                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
802                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
803                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
804                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
805
806             ENDDO
807
808!--          Linear interpolate to zw grid
809             DO k = nzb+2, nzt+8
810                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
811                                   rrtm_tlay(0,k-1))                           &
812                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
813                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
814             ENDDO
815
816
817!
818!--          Calculate liquid water path and cloud fraction for each column.
819!--          Note that LWP is required in g/m² instead of kg/kg m.
820             rrtm_cldfr  = 0.0_wp
821             rrtm_reliq  = 0.0_wp
822             rrtm_cliqwp = 0.0_wp
823             rrtm_icld   = 0
824
825             DO k = nzb+1, nzt+1
826                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                    &
827                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))        &
828                                    * 100.0_wp / g
829
830                IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
831                   rrtm_cldfr(0,k) = 1.0_wp
832                   IF ( rrtm_icld == 0 )  rrtm_icld = 1
833
834!
835!--                Calculate cloud droplet effective radius
836                   IF ( cloud_physics )  THEN
837                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
838                                        * rho_surface                          &
839                                        / ( 4.0_wp * pi * nc_const * rho_l )   &
840                                        )**0.33333333333333_wp                 &
841                                        * EXP( LOG( sigma_gc )**2 )
842
843                   ELSEIF ( cloud_droplets )  THEN
844                      number_of_particles = prt_count(k,j,i)
845
846                      IF (number_of_particles <= 0)  CYCLE
847                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
848                      s_r2 = 0.0_wp
849                      s_r3 = 0.0_wp
850
851                      DO  n = 1, number_of_particles
852                         IF ( particles(n)%particle_mask )  THEN
853                            s_r2 = s_r2 + particles(n)%radius**2 * &
854                                   particles(n)%weight_factor
855                            s_r3 = s_r3 + particles(n)%radius**3 * &
856                                   particles(n)%weight_factor
857                         ENDIF
858                      ENDDO
859
860                      IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
861
862                   ENDIF
863
864!
865!--                Limit effective radius
866                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
867                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
868                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
869                  ENDIF
870                ENDIF
871             ENDDO
872
873!
874!--          Set surface temperature
875             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
876
877             IF ( lw_radiation )  THEN
878               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
879               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
880               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
881               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
882               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
883               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
884               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
885               rrtm_reliq      , rrtm_lw_tauaer,                               &
886               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
887               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
888               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
889
890!
891!--             Save fluxes
892                DO k = nzb, nzt+1
893                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
894                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
895                ENDDO
896
897!
898!--             Save heating rates (convert from K/d to K/h)
899                DO k = nzb+1, nzt+1
900                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
901                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
902                ENDDO
903
904             ENDIF
905
906             IF ( sw_radiation .AND. sun_up )  THEN
907                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
908               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
909               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
910               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
911               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
912               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
913               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
914               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
915               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
916               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
917               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
918               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
919 
920!
921!--             Save fluxes
922                DO k = nzb, nzt+1
923                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
924                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
925                ENDDO
926
927!
928!--             Save heating rates (convert from K/d to K/s)
929                DO k = nzb+1, nzt+1
930                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
931                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
932                ENDDO
933
934             ENDIF
935
936!
937!--          Calculate surface net radiation
938             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
939                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
940
941          ENDDO
942       ENDDO
943
944       CALL exchange_horiz( rad_lw_in,  nbgp )
945       CALL exchange_horiz( rad_lw_out, nbgp )
946       CALL exchange_horiz( rad_lw_hr,    nbgp )
947       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
948
949       CALL exchange_horiz( rad_sw_in,  nbgp )
950       CALL exchange_horiz( rad_sw_out, nbgp ) 
951       CALL exchange_horiz( rad_sw_hr,    nbgp )
952       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
953
954       CALL exchange_horiz_2d( rad_net, nbgp )
955#endif
956
957    END SUBROUTINE radiation_rrtmg
958
959
960!------------------------------------------------------------------------------!
961! Description:
962! ------------
963!> Calculate the cosine of the zenith angle (variable is called zenith)
964!------------------------------------------------------------------------------!
965    SUBROUTINE calc_zenith
966
967       IMPLICIT NONE
968
969       REAL(wp) ::  declination,  & !< solar declination angle
970                    hour_angle      !< solar hour angle
971!
972!--    Calculate current day and time based on the initial values and simulation
973!--    time
974       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
975                               / 86400.0_wp ), KIND=iwp)
976       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
977
978
979!
980!--    Calculate solar declination and hour angle   
981       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
982       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
983
984!
985!--    Calculate zenith angle
986       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
987                                            * COS(hour_angle)
988       zenith(0) = MAX(0.0_wp,zenith(0))
989
990!
991!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
992       IF ( zenith(0) > 0.0_wp )  THEN
993          sun_up = .TRUE.
994       ELSE
995          sun_up = .FALSE.
996       END IF
997
998    END SUBROUTINE calc_zenith
999
1000#if defined ( __rrtmg ) && defined ( __netcdf )
1001!------------------------------------------------------------------------------!
1002! Description:
1003! ------------
1004!> Calculates surface albedo components based on Briegleb (1992) and
1005!> Briegleb et al. (1986)
1006!------------------------------------------------------------------------------!
1007    SUBROUTINE calc_albedo
1008
1009        IMPLICIT NONE
1010
1011        IF ( sun_up )  THEN
1012!
1013!--        Ocean
1014           IF ( albedo_type == 1 )  THEN
1015              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
1016                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
1017                                            * ( zenith(0) - 0.5_wp )           &
1018                                            * ( zenith(0) - 1.0_wp )
1019              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1020!
1021!--        Snow
1022           ELSEIF ( albedo_type == 16 )  THEN
1023              IF ( zenith(0) < 0.5_wp )  THEN
1024                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
1025                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1026                                     * zenith(0))) - 1.0_wp
1027                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
1028                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1029                                     * zenith(0))) - 1.0_wp
1030
1031                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1032                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1033              ELSE
1034                 rrtm_aldir(0,:,:) = aldif
1035                 rrtm_asdir(0,:,:) = asdif
1036              ENDIF
1037!
1038!--        Sea ice
1039           ELSEIF ( albedo_type == 15 )  THEN
1040                 rrtm_aldir(0,:,:) = aldif
1041                 rrtm_asdir(0,:,:) = asdif
1042!
1043!--        Land surfaces
1044           ELSE
1045              SELECT CASE ( albedo_type )
1046
1047!
1048!--              Surface types with strong zenith dependence
1049                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
1050                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
1051                                        (1.0_wp + 0.8_wp * zenith(0))
1052                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
1053                                        (1.0_wp + 0.8_wp * zenith(0))
1054!
1055!--              Surface types with weak zenith dependence
1056                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
1057                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
1058                                        (1.0_wp + 0.2_wp * zenith(0))
1059                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
1060                                        (1.0_wp + 0.2_wp * zenith(0))
1061
1062                 CASE DEFAULT
1063
1064              END SELECT
1065           ENDIF
1066!
1067!--        Diffusive albedo is taken from Table 2
1068           rrtm_aldif(0,:,:) = aldif
1069           rrtm_asdif(0,:,:) = asdif
1070
1071        ELSE
1072
1073           rrtm_aldir(0,:,:) = 0.0_wp
1074           rrtm_asdir(0,:,:) = 0.0_wp
1075           rrtm_aldif(0,:,:) = 0.0_wp
1076           rrtm_asdif(0,:,:) = 0.0_wp
1077        ENDIF
1078    END SUBROUTINE calc_albedo
1079
1080!------------------------------------------------------------------------------!
1081! Description:
1082! ------------
1083!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1084!------------------------------------------------------------------------------!
1085    SUBROUTINE read_sounding_data
1086
1087       USE netcdf_control
1088
1089       IMPLICIT NONE
1090
1091       INTEGER(iwp) :: id,           & !< NetCDF id of input file
1092                       id_dim_zrad,  & !< pressure level id in the NetCDF file
1093                       id_var,       & !< NetCDF variable id
1094                       k,            & !< loop index
1095                       nz_snd,       & !< number of vertical levels in the sounding data
1096                       nz_snd_start, & !< start vertical index for sounding data to be used
1097                       nz_snd_end      !< end vertical index for souding data to be used
1098
1099       REAL(wp) :: t_surface           !< actual surface temperature
1100
1101       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1102                                               t_snd_tmp      !< temporary temperature profile (sounding)
1103
1104!
1105!--    In case of updates, deallocate arrays first (sufficient to check one
1106!--    array as the others are automatically allocated). This is required
1107!--    because nzt_rad might change during the update
1108       IF ( ALLOCATED ( hyp_snd ) )  THEN
1109          DEALLOCATE( hyp_snd )
1110          DEALLOCATE( t_snd )
1111          DEALLOCATE( q_snd  )
1112          DEALLOCATE ( rrtm_play )
1113          DEALLOCATE ( rrtm_plev )
1114          DEALLOCATE ( rrtm_tlay )
1115          DEALLOCATE ( rrtm_tlev )
1116
1117          DEALLOCATE ( rrtm_h2ovmr )
1118          DEALLOCATE ( rrtm_cicewp )
1119          DEALLOCATE ( rrtm_cldfr )
1120          DEALLOCATE ( rrtm_cliqwp )
1121          DEALLOCATE ( rrtm_reice )
1122          DEALLOCATE ( rrtm_reliq )
1123          DEALLOCATE ( rrtm_lw_taucld )
1124          DEALLOCATE ( rrtm_lw_tauaer )
1125
1126          DEALLOCATE ( rrtm_lwdflx  )
1127          DEALLOCATE ( rrtm_lwdflxc )
1128          DEALLOCATE ( rrtm_lwuflx  )
1129          DEALLOCATE ( rrtm_lwuflxc )
1130          DEALLOCATE ( rrtm_lwuflx_dt )
1131          DEALLOCATE ( rrtm_lwuflxc_dt )
1132          DEALLOCATE ( rrtm_lwhr  )
1133          DEALLOCATE ( rrtm_lwhrc )
1134
1135          DEALLOCATE ( rrtm_sw_taucld )
1136          DEALLOCATE ( rrtm_sw_ssacld )
1137          DEALLOCATE ( rrtm_sw_asmcld )
1138          DEALLOCATE ( rrtm_sw_fsfcld )
1139          DEALLOCATE ( rrtm_sw_tauaer )
1140          DEALLOCATE ( rrtm_sw_ssaaer )
1141          DEALLOCATE ( rrtm_sw_asmaer ) 
1142          DEALLOCATE ( rrtm_sw_ecaer )   
1143 
1144          DEALLOCATE ( rrtm_swdflx  )
1145          DEALLOCATE ( rrtm_swdflxc )
1146          DEALLOCATE ( rrtm_swuflx  )
1147          DEALLOCATE ( rrtm_swuflxc )
1148          DEALLOCATE ( rrtm_swhr  )
1149          DEALLOCATE ( rrtm_swhrc )
1150
1151       ENDIF
1152
1153!
1154!--    Open file for reading
1155       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1156       CALL handle_netcdf_error( 'netcdf', 549 )
1157
1158!
1159!--    Inquire dimension of z axis and save in nz_snd
1160       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1161       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1162       CALL handle_netcdf_error( 'netcdf', 551 )
1163
1164!
1165! !--    Allocate temporary array for storing pressure data
1166       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1167       hyp_snd_tmp = 0.0_wp
1168
1169
1170!--    Read pressure from file
1171       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1172       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
1173                               count = (/nz_snd/) )
1174       CALL handle_netcdf_error( 'netcdf', 552 )
1175
1176!
1177!--    Allocate temporary array for storing temperature data
1178       ALLOCATE( t_snd_tmp(1:nz_snd) )
1179       t_snd_tmp = 0.0_wp
1180
1181!
1182!--    Read temperature from file
1183       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1184       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
1185                               count = (/nz_snd/) )
1186       CALL handle_netcdf_error( 'netcdf', 553 )
1187
1188!
1189!--    Calculate start of sounding data
1190       nz_snd_start = nz_snd + 1
1191       nz_snd_end   = nz_snd + 1
1192
1193!
1194!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
1195!--    in Pa, hyp_snd in hPa).
1196       DO  k = 1, nz_snd
1197          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
1198             nz_snd_start = k
1199             EXIT
1200          END IF
1201       END DO
1202
1203       IF ( nz_snd_start <= nz_snd )  THEN
1204          nz_snd_end = nz_snd
1205       END IF
1206
1207
1208!
1209!--    Calculate of total grid points for RRTMG calculations
1210       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1211
1212!
1213!--    Save data above LES domain in hyp_snd, t_snd and q_snd
1214!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1215       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1216       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
1217       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
1218       hyp_snd = 0.0_wp
1219       t_snd = 0.0_wp
1220       q_snd = 0.0_wp
1221
1222       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start:nz_snd_end)
1223       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start:nz_snd_end)
1224
1225       DEALLOCATE ( hyp_snd_tmp )
1226       DEALLOCATE ( t_snd_tmp )
1227
1228       nc_stat = NF90_CLOSE( id )
1229
1230!
1231!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
1232!--    top of the LES domain. This routine does not consider horizontal or
1233!--    vertical variability of pressure and temperature
1234       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
1235       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
1236
1237       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1238       DO k = nzb+1, nzt+1
1239          rrtm_play(0,k) = hyp(k) * 0.01_wp
1240          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1241                         t_surface )**(1.0_wp/0.286_wp)
1242       ENDDO
1243
1244       DO k = nzt+2, nzt_rad
1245          rrtm_play(0,k) = hyp_snd(k)
1246          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1247       ENDDO
1248       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
1249                                   1.5 * hyp_snd(nzt_rad)                      &
1250                                 - 0.5 * hyp_snd(nzt_rad-1) )
1251       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
1252                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
1253
1254       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1255
1256!
1257!--    Calculate temperature/humidity levels at top of the LES domain.
1258!--    Currently, the temperature is taken from sounding data (might lead to a
1259!--    temperature jump at interface. To do: Humidity is currently not
1260!--    calculated above the LES domain.
1261       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
1262       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
1263       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1264
1265       DO k = nzt+8, nzt_rad
1266          rrtm_tlay(0,k)   = t_snd(k)
1267          rrtm_h2ovmr(0,k) = q_snd(k)
1268       ENDDO
1269       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
1270                                - rrtm_tlay(0,nzt_rad-1)
1271       DO k = nzt+9, nzt_rad+1
1272          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
1273                             - rrtm_tlay(0,k-1))                               &
1274                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
1275                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1276       ENDDO
1277       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1278
1279       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
1280                                  - rrtm_tlev(0,nzt_rad)
1281!
1282!--    Allocate remaining RRTMG arrays
1283       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1284       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1285       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1286       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1287       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1288       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1289       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1290       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1291       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1292       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1293       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1294       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1295       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1296       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
1297       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
1298
1299!
1300!--    The ice phase is currently not considered in PALM
1301       rrtm_cicewp = 0.0_wp
1302       rrtm_reice  = 0.0_wp
1303
1304!
1305!--    Set other parameters (move to NAMELIST parameters in the future)
1306       rrtm_lw_tauaer = 0.0_wp
1307       rrtm_lw_taucld = 0.0_wp
1308       rrtm_sw_taucld = 0.0_wp
1309       rrtm_sw_ssacld = 0.0_wp
1310       rrtm_sw_asmcld = 0.0_wp
1311       rrtm_sw_fsfcld = 0.0_wp
1312       rrtm_sw_tauaer = 0.0_wp
1313       rrtm_sw_ssaaer = 0.0_wp
1314       rrtm_sw_asmaer = 0.0_wp
1315       rrtm_sw_ecaer  = 0.0_wp
1316
1317
1318       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
1319       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
1320       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
1321       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1322       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1323       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1324
1325       rrtm_swdflx  = 0.0_wp
1326       rrtm_swuflx  = 0.0_wp
1327       rrtm_swhr    = 0.0_wp 
1328       rrtm_swuflxc = 0.0_wp
1329       rrtm_swdflxc = 0.0_wp
1330       rrtm_swhrc   = 0.0_wp
1331
1332       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
1333       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
1334       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
1335       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1336       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1337       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1338
1339       rrtm_lwdflx  = 0.0_wp
1340       rrtm_lwuflx  = 0.0_wp
1341       rrtm_lwhr    = 0.0_wp 
1342       rrtm_lwuflxc = 0.0_wp
1343       rrtm_lwdflxc = 0.0_wp
1344       rrtm_lwhrc   = 0.0_wp
1345
1346       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
1347       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
1348
1349       rrtm_lwuflxc_dt = 0.0_wp
1350       rrtm_lwuflxc_dt = 0.0_wp
1351
1352    END SUBROUTINE read_sounding_data
1353
1354
1355!------------------------------------------------------------------------------!
1356! Description:
1357! ------------
1358!> Read trace gas data from file
1359!------------------------------------------------------------------------------!
1360    SUBROUTINE read_trace_gas_data
1361
1362       USE netcdf_control
1363       USE rrsw_ncpar
1364
1365       IMPLICIT NONE
1366
1367       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
1368
1369       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
1370           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
1371                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
1372
1373       INTEGER(iwp) :: id,     & !< NetCDF id
1374                       k,      & !< loop index
1375                       m,      & !< loop index
1376                       n,      & !< loop index
1377                       nabs,   & !< number of absorbers
1378                       np,     & !< number of pressure levels
1379                       id_abs, & !< NetCDF id of the respective absorber
1380                       id_dim, & !< NetCDF id of asborber's dimension
1381                       id_var    !< NetCDf id ot the absorber
1382
1383       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
1384
1385
1386       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
1387                                                 rrtm_play_tmp, & !< temporary array for pressure zu-levels
1388                                                 rrtm_plev_tmp, & !< temporary array for pressure zw-levels
1389                                                 trace_path_tmp   !< temporary array for storing trace gas path data
1390
1391       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
1392                                                 trace_mls_path, & !< array for storing trace gas path data
1393                                                 trace_mls_tmp     !< temporary array for storing trace gas data
1394
1395
1396!
1397!--    In case of updates, deallocate arrays first (sufficient to check one
1398!--    array as the others are automatically allocated)
1399       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
1400          DEALLOCATE ( rrtm_o3vmr  )
1401          DEALLOCATE ( rrtm_co2vmr )
1402          DEALLOCATE ( rrtm_ch4vmr )
1403          DEALLOCATE ( rrtm_n2ovmr )
1404          DEALLOCATE ( rrtm_o2vmr  )
1405          DEALLOCATE ( rrtm_cfc11vmr )
1406          DEALLOCATE ( rrtm_cfc12vmr )
1407          DEALLOCATE ( rrtm_cfc22vmr )
1408          DEALLOCATE ( rrtm_ccl4vmr  )
1409       ENDIF
1410
1411!
1412!--    Allocate trace gas profiles
1413       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
1414       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
1415       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
1416       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
1417       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
1418       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
1419       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
1420       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
1421       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
1422
1423!
1424!--    Open file for reading
1425       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1426       CALL handle_netcdf_error( 'netcdf', 549 )
1427!
1428!--    Inquire dimension ids and dimensions
1429       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
1430       CALL handle_netcdf_error( 'netcdf', 550 )
1431       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
1432       CALL handle_netcdf_error( 'netcdf', 550 )
1433
1434       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
1435       CALL handle_netcdf_error( 'netcdf', 550 )
1436       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
1437       CALL handle_netcdf_error( 'netcdf', 550 )
1438   
1439
1440!
1441!--    Allocate pressure, and trace gas arrays     
1442       ALLOCATE( p_mls(1:np) )
1443       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
1444       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
1445
1446
1447       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1448       CALL handle_netcdf_error( 'netcdf', 550 )
1449       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
1450       CALL handle_netcdf_error( 'netcdf', 550 )
1451
1452       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
1453       CALL handle_netcdf_error( 'netcdf', 550 )
1454       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
1455       CALL handle_netcdf_error( 'netcdf', 550 )
1456
1457
1458!
1459!--    Write absorber amounts (mls) to trace_mls
1460       DO n = 1, num_trace_gases
1461          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
1462
1463          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
1464
1465!
1466!--       Replace missing values by zero
1467          WHERE ( trace_mls(n,:) > 2.0_wp ) 
1468             trace_mls(n,:) = 0.0_wp
1469          END WHERE
1470       END DO
1471
1472       DEALLOCATE ( trace_mls_tmp )
1473
1474       nc_stat = NF90_CLOSE( id )
1475       CALL handle_netcdf_error( 'netcdf', 551 )
1476
1477!
1478!--    Add extra pressure level for calculations of the trace gas paths
1479       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
1480       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
1481
1482       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
1483       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
1484       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
1485       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
1486                                         * rrtm_plev(0,nzt_rad+1) )
1487 
1488!
1489!--    Calculate trace gas path (zero at surface) with interpolation to the
1490!--    sounding levels
1491       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
1492
1493       trace_mls_path(nzb+1,:) = 0.0_wp
1494       
1495       DO k = nzb+2, nzt_rad+2
1496          DO m = 1, num_trace_gases
1497             trace_mls_path(k,m) = trace_mls_path(k-1,m)
1498
1499!
1500!--          When the pressure level is higher than the trace gas pressure
1501!--          level, assume that
1502             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
1503               
1504                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
1505                                      * ( rrtm_plev_tmp(k-1)                   &
1506                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
1507                                        ) / g
1508             ENDIF
1509
1510!
1511!--          Integrate for each sounding level from the contributing p_mls
1512!--          levels
1513             DO n = 2, np
1514!
1515!--             Limit p_mls so that it is within the model level
1516                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
1517                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
1518                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
1519                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
1520
1521                IF ( p_mls_l > p_mls_u )  THEN
1522
1523!
1524!--                Calculate weights for interpolation
1525                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
1526                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
1527                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
1528
1529!
1530!--                Add level to trace gas path
1531                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
1532                                         +  ( p_wgt_u * trace_mls(m,n)         &
1533                                            + p_wgt_l * trace_mls(m,n-1) )     &
1534                                         * (p_mls_l - p_mls_u) / g
1535                ENDIF
1536             ENDDO
1537
1538             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
1539                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
1540                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
1541                                          - rrtm_plev_tmp(k)                   &
1542                                        ) / g
1543             ENDIF 
1544          ENDDO
1545       ENDDO
1546
1547
1548!
1549!--    Prepare trace gas path profiles
1550       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
1551
1552       DO m = 1, num_trace_gases
1553
1554          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
1555                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
1556                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
1557                                       - rrtm_plev_tmp(2:nzt_rad+2) )
1558
1559!
1560!--       Save trace gas paths to the respective arrays
1561          SELECT CASE ( TRIM( trace_names(m) ) )
1562
1563             CASE ( 'O3' )
1564
1565                rrtm_o3vmr(0,:) = trace_path_tmp(:)
1566
1567             CASE ( 'CO2' )
1568
1569                rrtm_co2vmr(0,:) = trace_path_tmp(:)
1570
1571             CASE ( 'CH4' )
1572
1573                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
1574
1575             CASE ( 'N2O' )
1576
1577                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
1578
1579             CASE ( 'O2' )
1580
1581                rrtm_o2vmr(0,:) = trace_path_tmp(:)
1582
1583             CASE ( 'CFC11' )
1584
1585                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
1586
1587             CASE ( 'CFC12' )
1588
1589                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
1590
1591             CASE ( 'CFC22' )
1592
1593                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
1594
1595             CASE ( 'CCL4' )
1596
1597                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
1598
1599             CASE DEFAULT
1600
1601          END SELECT
1602
1603       ENDDO
1604
1605       DEALLOCATE ( trace_path_tmp )
1606       DEALLOCATE ( trace_mls_path )
1607       DEALLOCATE ( rrtm_play_tmp )
1608       DEALLOCATE ( rrtm_plev_tmp )
1609       DEALLOCATE ( trace_mls )
1610       DEALLOCATE ( p_mls )
1611
1612    END SUBROUTINE read_trace_gas_data
1613
1614#endif
1615
1616
1617!------------------------------------------------------------------------------!
1618! Description:
1619! ------------
1620!> Calculate temperature tendency due to radiative cooling/heating.
1621!> Cache-optimized version.
1622!------------------------------------------------------------------------------!
1623    SUBROUTINE radiation_tendency_ij ( i, j, tend )
1624
1625       USE cloud_parameters,                                                   &
1626           ONLY:  pt_d_t
1627
1628       IMPLICIT NONE
1629
1630       INTEGER(iwp) :: i, j, k !< loop indices
1631
1632       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
1633
1634#if defined ( __rrtmg )
1635!
1636!--    Calculate tendency based on heating rate
1637       DO k = nzb+1, nzt+1
1638          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
1639                                      * pt_d_t(k) * d_seconds_hour
1640       ENDDO
1641
1642#endif
1643
1644    END SUBROUTINE radiation_tendency_ij
1645
1646
1647!------------------------------------------------------------------------------!
1648! Description:
1649! ------------
1650!> Calculate temperature tendency due to radiative cooling/heating.
1651!> Vector-optimized version
1652!------------------------------------------------------------------------------!
1653    SUBROUTINE radiation_tendency ( tend )
1654
1655       USE cloud_parameters,                                                   &
1656           ONLY:  pt_d_t
1657
1658       USE indices,                                                            &
1659           ONLY:  nxl, nxr, nyn, nys
1660
1661       IMPLICIT NONE
1662
1663       INTEGER(iwp) :: i, j, k !< loop indices
1664
1665       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
1666
1667#if defined ( __rrtmg )
1668!
1669!--    Calculate tendency based on heating rate
1670       DO  i = nxl, nxr
1671          DO  j = nys, nyn
1672             DO k = nzb+1, nzt+1
1673                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
1674                                            +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
1675                                            * d_seconds_hour
1676             ENDDO
1677         ENDDO
1678       ENDDO
1679#endif
1680
1681    END SUBROUTINE radiation_tendency
1682
1683 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.