Changeset 1585 for palm/trunk/SOURCE/radiation_model.f90
- Timestamp:
- Apr 30, 2015 7:05:52 AM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model.f90
r1572 r1585 15 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 16 16 ! 17 ! Copyright 1997-201 4Leibniz Universitaet Hannover17 ! Copyright 1997-2015 Leibniz Universitaet Hannover 18 18 !--------------------------------------------------------------------------------! 19 19 ! 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Added support for RRTMG 23 23 ! 24 24 ! Former revisions: … … 40 40 ! Description: 41 41 ! ------------ 42 ! Radiation model(s), to be used e.g. with the land surface scheme 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) 43 53 !------------------------------------------------------------------------------! 44 54 45 55 USE arrays_3d, & 46 ONLY: pt 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 47 63 48 64 USE control_parameters, & 49 ONLY: phi, surface_pressure, time_since_reference_point 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 50 68 51 69 USE indices, & 52 ONLY: nxl g, nxrg, nyng, nysg, nzb_s_inner70 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt 53 71 54 72 USE kinds 73 74 USE netcdf 55 75 56 76 USE netcdf_control, & 57 77 ONLY: dots_label, dots_num, dots_unit 58 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 59 100 60 101 IMPLICIT NONE 61 102 62 CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm' 63 64 INTEGER(iwp) :: i, j, k 65 66 67 INTEGER(iwp) :: day_init = 172, & !: day of the year at model start (21/06) 68 dots_rad = 0, & !: starting index for timeseries output 69 irad_scheme = 0 70 71 LOGICAL :: radiation = .FALSE. !: flag parameter indicating wheather the radiation model is used 72 73 REAL(wp), PARAMETER :: sw_0 = 1368.0_wp, & !: solar constant 74 pi = 3.14159265358979323_wp, & 75 sigma_sb = 5.67E-8_wp !: Stefan-Boltzmann constant 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 76 145 77 REAL(wp) :: albedo = 0.2_wp, & !: NAMELIST alpha 78 dt_radiation = 0.0_wp, & !: radiation model timestep 79 exn, & !: Exner function 80 lon = 0.0_wp, & !: longitude in radians 81 lat = 0.0_wp, & !: latitude in radians 82 decl_1, & !: declination coef. 1 83 decl_2, & !: declination coef. 2 84 decl_3, & !: declination coef. 3 85 time_utc, & !: current time in UTC 86 time_utc_init = 43200.0_wp, & !: UTC time at model start (noon) 87 day, & !: current day of the year 88 lambda = 0.0_wp, & !: longitude in degrees 89 declination, & !: solar declination angle 90 net_radiation = 0.0_wp, & !: net radiation at surface 91 hour_angle, & !: solar hour angle 92 time_radiation = 0.0_wp, & 93 zenith, & !: solar zenith angle 94 sky_trans !: sky transmissivity 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 95 168 96 169 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 97 alpha, & !: surface albedo170 alpha, & !: surface broadband albedo (used for clear-sky scheme) 98 171 rad_net, & !: net radiation at the surface 99 rad_net_av, & !: average of rad_net 100 rad_lw_in, & !: incoming longwave radiation 101 rad_lw_out, & !: outgoing longwave radiation 102 rad_sw_in, & !: incoming shortwave radiation 103 rad_sw_in_av, & !: average of rad_sw_in 104 rad_sw_out !: outgoing shortwave radiation 105 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 106 300 107 301 INTERFACE init_radiation … … 113 307 END INTERFACE radiation_clearsky 114 308 115 INTERFACE radiation_ constant116 MODULE PROCEDURE radiation_ constant117 END INTERFACE radiation_ constant118 119 INTERFACE radiation_ rrtm120 MODULE PROCEDURE radiation_ rrtm121 END INTERFACE radiation_rrtm122 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 123 317 124 318 SAVE … … 126 320 PRIVATE 127 321 128 PUBLIC albedo, day_init, dots_rad, dt_radiation, init_radiation, & 129 irad_scheme, lambda, net_radiation, rad_net, rad_net_av, radiation, & 130 radiation_clearsky, radiation_constant, radiation_rrtm, & 131 radiation_scheme, rad_sw_in, rad_sw_in_av, sigma_sb, & 132 time_radiation, time_utc_init 133 134 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 135 334 136 335 CONTAINS … … 143 342 SUBROUTINE init_radiation 144 343 145 146 344 IMPLICIT NONE 147 345 148 ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) ) 149 ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) ) 150 ALLOCATE ( rad_lw_in(nysg:nyng,nxlg:nxrg) ) 151 ALLOCATE ( rad_lw_out(nysg:nyng,nxlg:nxrg) ) 152 ALLOCATE ( rad_sw_in(nysg:nyng,nxlg:nxrg) ) 153 ALLOCATE ( rad_sw_out(nysg:nyng,nxlg:nxrg) ) 154 155 rad_sw_in = 0.0_wp 156 rad_sw_out = 0.0_wp 157 rad_lw_in = 0.0_wp 158 rad_lw_out = 0.0_wp 159 rad_net = 0.0_wp 160 161 alpha = albedo 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 162 352 163 353 ! 164 354 !-- Fix net radiation in case of radiation_scheme = 'constant' 165 IF ( irad_scheme == 0) THEN355 IF ( radiation_scheme == 'constant' ) THEN 166 356 rad_net = net_radiation 167 ! 168 !-- Calculate radiation scheme constants 169 ELSEIF ( irad_scheme == 1 .OR. irad_scheme == 2 ) THEN 357 radiation = .FALSE. 358 ! 359 !-- Calculate orbital constants 360 ELSE 170 361 decl_1 = SIN(23.45_wp * pi / 180.0_wp) 171 362 decl_2 = 2.0_wp * pi / 365.0_wp 172 363 decl_3 = decl_2 * 81.0_wp 173 ! 174 !-- Calculate latitude and longitude angles (lat and lon, respectively) 175 lat = phi * pi / 180.0_wp 176 lon = lambda * pi / 180.0_wp 364 lat = phi * pi / 180.0_wp 365 lon = lambda * pi / 180.0_wp 177 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 178 558 ! 179 559 !-- Add timeseries for radiation model 560 dots_rad = dots_num + 1 180 561 dots_label(dots_num+1) = "rad_net" 181 dots_label(dots_num+2) = "rad_sw_in" 182 dots_unit(dots_num+1:dots_num+2) = "W/m2" 183 184 dots_rad = dots_num + 1 185 dots_num = dots_num + 2 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 186 589 187 590 RETURN … … 196 599 !------------------------------------------------------------------------------! 197 600 SUBROUTINE radiation_clearsky 198 601 602 USE indices, & 603 ONLY: nbgp 199 604 200 605 IMPLICIT NONE 201 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 202 859 ! 203 860 !-- Calculate current day and time based on the initial values and simulation 204 861 !-- time 205 day = day_init + FLOOR( (time_utc_init + time_since_reference_point) &206 / 86400.0_wp ) 862 day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) & 863 / 86400.0_wp ), KIND=iwp) 207 864 time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp) 208 865 … … 210 867 ! 211 868 !-- Calculate solar declination and hour angle 212 declination = ASIN( decl_1 * SIN(decl_2 * day- decl_3) )869 declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) ) 213 870 hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi 214 871 215 872 ! 216 873 !-- Calculate zenith angle 217 zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)&874 zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) & 218 875 * COS(hour_angle) 219 zenith = MAX(0.0_wp,zenith) 220 221 ! 222 !-- Calculate sky transmissivity 223 sky_trans = 0.6_wp + 0.2_wp * zenith 224 225 ! 226 !-- Calculate value of the Exner function 227 exn = (surface_pressure / 1000.0_wp )**0.286_wp 228 229 ! 230 !-- Calculate radiation fluxes and net radiation (rad_net) for each grid 231 !-- point 232 DO i = nxlg, nxrg 233 DO j = nysg, nyng 234 235 k = nzb_s_inner(j,i) 236 rad_sw_in(j,i) = sw_0 * sky_trans * zenith 237 rad_sw_out(j,i) = - alpha(j,i) * rad_sw_in(j,i) 238 rad_lw_out(j,i) = - sigma_sb * (pt(k,j,i) * exn)**4 239 rad_lw_in(j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4 240 rad_net(j,i) = rad_sw_in(j,i) + rad_sw_out(j,i) & 241 + rad_lw_in(j,i) + rad_lw_out(j,i) 242 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 243 1405 ENDDO 244 1406 ENDDO 245 1407 246 RETURN 247 248 END SUBROUTINE radiation_clearsky 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 249 1477 250 1478 !------------------------------------------------------------------------------! 251 1479 ! Description: 252 1480 ! ------------ 253 !-- Prescribed net radiation 254 !------------------------------------------------------------------------------! 255 SUBROUTINE radiation_constant 256 257 rad_net = net_radiation 258 259 END SUBROUTINE radiation_constant 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 260 1528 261 1529 !------------------------------------------------------------------------------! 262 1530 ! Description: 263 1531 ! ------------ 264 !-- Implementation of the RRTM radiation_scheme 265 !------------------------------------------------------------------------------! 266 SUBROUTINE radiation_rrtm 267 268 269 END SUBROUTINE radiation_rrtm 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(k ,j,i) - rad_sw_out(k ,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(k ,j,i) - rad_lw_out(k ,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 270 1585 271 1586 END MODULE radiation_model_mod
Note: See TracChangeset
for help on using the changeset viewer.
