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

Last change on this file since 1590 was 1590, checked in by maronga, 7 years ago

bugfix in definition of character arrays within the LSM and radiation code

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