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

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

bugfix: radiation_scheme = constant caused model crash when used in combination with land surface model

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