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

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

Added support for RRTMG radiation code

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