Changeset 3824 for palm/trunk/SOURCE
- Timestamp:
- Mar 27, 2019 3:56:16 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
chem_photolysis_mod.f90 (modified) (3 diffs)
-
chemistry_model_mod.f90 (modified) (7 diffs)
-
radiation_model_mod.f90 (modified) (30 diffs)
-
urban_surface_mod.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/chem_photolysis_mod.f90
r3614 r3824 234 234 235 235 USE radiation_model_mod, & 236 ONLY: calc_zenith, zenith236 ONLY: calc_zenith, cos_zenith 237 237 238 238 IMPLICIT NONE … … 247 247 CALL calc_zenith 248 248 249 IF ( zenith(0)> 0.0_wp ) THEN250 coszi = 1. / zenith(0)249 IF ( cos_zenith > 0.0_wp ) THEN 250 coszi = 1. / cos_zenith 251 251 252 252 DO iphot = 1,nphot … … 254 254 IF ( TRIM( names_s(iav) ) == TRIM( phot_names(iphot) ) ) then 255 255 phot_frequen(iphot)%freq(nzb+1:nzt,:,:) = & 256 par_l(iav) * zenith(0)**par_m(iav) * EXP( -par_n(iav) * coszi )256 par_l(iav) * cos_zenith**par_m(iav) * EXP( -par_n(iav) * coszi ) 257 257 ENDIF 258 258 ENDDO -
palm/trunk/SOURCE/chemistry_model_mod.f90
r3821 r3824 2799 2799 2800 2800 USE radiation_model_mod, & 2801 ONLY: zenith2801 ONLY: cos_zenith 2802 2802 2803 2803 … … 3419 3419 ! 3420 3420 !-- Get Rc_tot 3421 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &3421 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3422 3422 relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3423 3423 r_aero_lu , Rb ) … … 3566 3566 ! 3567 3567 !-- Get Rc_tot 3568 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &3568 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3569 3569 relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3570 3570 r_aero_lu , Rb ) … … 3713 3713 ! 3714 3714 !-- Get Rc_tot 3715 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &3715 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 3716 3716 relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 3717 3717 r_aero_lu , Rb ) … … 4109 4109 ! 4110 4110 !-- Get Rc_tot 4111 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &4111 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4112 4112 relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4113 4113 r_aero_lu , Rb ) … … 4249 4249 ! 4250 4250 !-- Get Rc_tot 4251 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &4251 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4252 4252 relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4253 4253 r_aero_lu , Rb ) … … 4395 4395 ! 4396 4396 !-- Get Rc_tot 4397 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &4397 CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, & 4398 4398 relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, & 4399 4399 r_aero_lu , Rb ) -
palm/trunk/SOURCE/radiation_model_mod.f90
r3814 r3824 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Change zenith(0:0) and others to scalar. 31 ! Code review. 32 ! 33 ! 3814 2019-03-26 08:40:31Z pavelkrc 30 34 ! Rename exported nzu, nzp and related variables due to name conflict 31 35 ! … … 747 751 748 752 749 REAL(wp) , DIMENSION(0:0) :: zenith, & !< cosine of solar zenith angle750 sun_dir_lat, & !< solar directional vector in latitudes751 sun_dir_lon !< solar directional vector in longitudes753 REAL(wp) :: cos_zenith !< cosine of solar zenith angle, also z-coordinate of solar unit vector 754 REAL(wp) :: sun_dir_lat !< y-coordinate of solar unit vector 755 REAL(wp) :: sun_dir_lon !< x-coordinate of solar unit vector 752 756 753 757 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_net_av !< average of net radiation (rad_net) at surface … … 1294 1298 rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, & 1295 1299 rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, & 1296 rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, s igma_sb, solar_constant,&1300 rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant, & 1297 1301 skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,& 1298 zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,&1302 cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, & 1299 1303 idir, jdir, kdir, id, iz, iy, ix, & 1300 1304 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & … … 2834 2838 ! 2835 2839 !-- Calculate sky transmissivity 2836 sky_trans = 0.6_wp + 0.2_wp * zenith(0)2840 sky_trans = 0.6_wp + 0.2_wp * cos_zenith 2837 2841 2838 2842 ! … … 2909 2913 k = nz_urban_t 2910 2914 2911 surf%rad_sw_in = solar_constant * sky_trans * zenith(0)2915 surf%rad_sw_in = solar_constant * sky_trans * cos_zenith 2912 2916 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 2913 2917 … … 2933 2937 k = surf%k(m) 2934 2938 2935 surf%rad_sw_in(m) = solar_constant * sky_trans * zenith(0)2939 surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith 2936 2940 2937 2941 ! … … 3379 3383 3380 3384 REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av 3385 REAL(wp), DIMENSION(0:0) :: zenith !< to provide indexed array 3381 3386 ! 3382 3387 !-- Just dummy arguments … … 3395 3400 !-- Calculate current (cosine of) zenith angle and whether the sun is up 3396 3401 CALL calc_zenith 3402 zenith(0) = cos_zenith 3397 3403 ! 3398 3404 !-- Calculate surface albedo. In case average radiation is applied, … … 4172 4178 ! 4173 4179 !-- Calculate cosine of solar zenith angle 4174 zenith(0)= SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &4180 cos_zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) & 4175 4181 * COS(hour_angle) 4176 zenith(0) = MAX(0.0_wp,zenith(0))4182 cos_zenith = MAX(0.0_wp,cos_zenith) 4177 4183 4178 4184 ! … … 4182 4188 ! 4183 4189 !-- Direction in longitudes equals to sin(solar_azimuth) * sin(zenith) 4184 sun_dir_lon (0)= -SIN(hour_angle) * COS(declination)4190 sun_dir_lon = -SIN(hour_angle) * COS(declination) 4185 4191 4186 4192 ! 4187 4193 !-- Direction in latitues equals to cos(solar_azimuth) * sin(zenith) 4188 sun_dir_lat (0)= SIN(declination) * COS(lat) - COS(hour_angle) &4194 sun_dir_lat = SIN(declination) * COS(lat) - COS(hour_angle) & 4189 4195 * COS(declination) * SIN(lat) 4190 4196 ENDIF … … 4192 4198 ! 4193 4199 !-- Check if the sun is up (otheriwse shortwave calculations can be skipped) 4194 IF ( zenith(0)> 0.0_wp ) THEN4200 IF ( cos_zenith > 0.0_wp ) THEN 4195 4201 sun_up = .TRUE. 4196 4202 ELSE … … 4227 4233 IF ( surf%albedo_type(ind_type,m) == 1 ) THEN 4228 4234 surf%rrtm_aldir(ind_type,m) = 0.026_wp / & 4229 ( zenith(0)**1.7_wp + 0.065_wp )&4230 + 0.15_wp * ( zenith(0)- 0.1_wp ) &4231 * ( zenith(0)- 0.5_wp ) &4232 * ( zenith(0)- 1.0_wp )4235 ( cos_zenith**1.7_wp + 0.065_wp )& 4236 + 0.15_wp * ( cos_zenith - 0.1_wp ) & 4237 * ( cos_zenith - 0.5_wp ) & 4238 * ( cos_zenith - 1.0_wp ) 4233 4239 surf%rrtm_asdir(ind_type,m) = surf%rrtm_aldir(ind_type,m) 4234 4240 ! 4235 4241 !-- Snow 4236 4242 ELSEIF ( surf%albedo_type(ind_type,m) == 16 ) THEN 4237 IF ( zenith(0)< 0.5_wp ) THEN4243 IF ( cos_zenith < 0.5_wp ) THEN 4238 4244 surf%rrtm_aldir(ind_type,m) = & 4239 4245 0.5_wp * ( 1.0_wp - surf%aldif(ind_type,m) ) & 4240 4246 * ( 3.0_wp / ( 1.0_wp + 4.0_wp & 4241 * zenith(0)) ) - 1.0_wp4247 * cos_zenith ) ) - 1.0_wp 4242 4248 surf%rrtm_asdir(ind_type,m) = & 4243 4249 0.5_wp * ( 1.0_wp - surf%asdif(ind_type,m) ) & 4244 4250 * ( 3.0_wp / ( 1.0_wp + 4.0_wp & 4245 * zenith(0)) ) - 1.0_wp4251 * cos_zenith ) ) - 1.0_wp 4246 4252 4247 4253 surf%rrtm_aldir(ind_type,m) = & … … 4282 4288 surf%rrtm_aldir(ind_type,m) = & 4283 4289 surf%aldif(ind_type,m) * 1.4_wp / & 4284 ( 1.0_wp + 0.8_wp * zenith(0))4290 ( 1.0_wp + 0.8_wp * cos_zenith ) 4285 4291 surf%rrtm_asdir(ind_type,m) = & 4286 4292 surf%asdif(ind_type,m) * 1.4_wp / & 4287 ( 1.0_wp + 0.8_wp * zenith(0))4293 ( 1.0_wp + 0.8_wp * cos_zenith ) 4288 4294 ! 4289 4295 !-- Surface types with weak zenith dependence … … 4291 4297 surf%rrtm_aldir(ind_type,m) = & 4292 4298 surf%aldif(ind_type,m) * 1.1_wp / & 4293 ( 1.0_wp + 0.2_wp * zenith(0))4299 ( 1.0_wp + 0.2_wp * cos_zenith ) 4294 4300 surf%rrtm_asdir(ind_type,m) = & 4295 4301 surf%asdif(ind_type,m) * 1.1_wp / & 4296 ( 1.0_wp + 0.2_wp * zenith(0))4302 ( 1.0_wp + 0.2_wp * cos_zenith ) 4297 4303 4298 4304 CASE DEFAULT … … 5018 5024 mrot(2, :) = (/ 0._wp, COS(alpha), SIN(alpha) /) 5019 5025 mrot(3, :) = (/ 0._wp, -SIN(alpha), COS(alpha) /) 5020 sunorig = (/ zenith(0), sun_dir_lat, sun_dir_lon /)5026 sunorig = (/ cos_zenith, sun_dir_lat, sun_dir_lon /) 5021 5027 sunorig = MATMUL(mrot, sunorig) 5022 5028 DO d = 0, nsurf_type … … 5024 5030 ENDDO 5025 5031 5026 IF ( zenith(0)> 0 ) THEN5032 IF ( cos_zenith > 0 ) THEN 5027 5033 !-- now we will "squash" the sunorig vector by grid box size in 5028 5034 !-- each dimension, so that this new direction vector will allow us … … 5186 5192 5187 5193 !-- direct radiation 5188 IF ( zenith(0)> 0 ) THEN5194 IF ( cos_zenith > 0 ) THEN 5189 5195 !--Identify solar direction vector (discretized number) 1) 5190 5196 !-- 5191 j = FLOOR(ACOS( zenith(0)) / pi * raytrace_discrete_elevs)5192 i = MODULO(NINT(ATAN2(sun_dir_lon (0), sun_dir_lat(0)) &5197 j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs) 5198 i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat) & 5193 5199 / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), & 5194 5200 raytrace_discrete_azims) … … 5199 5205 i = surfl(ix, isurf) 5200 5206 surfinswdir(isurf) = rad_sw_in_dir(j,i) * & 5201 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)5207 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / cos_zenith 5202 5208 ENDDO 5203 5209 ! … … 5207 5213 i = mrtbl(ix, imrt) 5208 5214 mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) & 5209 / zenith(0)/ 4._wp ! normal to sphere5215 / cos_zenith / 4._wp ! normal to sphere 5210 5216 ENDDO 5211 5217 ENDIF … … 5217 5223 mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc) 5218 5224 ENDDO 5219 5225 ! 5226 !-- Absorption in each local plant canopy grid box from the first atmospheric 5227 !-- pass of radiation 5220 5228 IF ( npcbl > 0 ) THEN 5221 5229 … … 5223 5231 pcbinswdif(:) = 0._wp 5224 5232 pcbinlw(:) = 0._wp 5225 ! 5226 !-- pcsf first pass 5233 5227 5234 DO icsf = 1, ncsfl 5228 5235 ipcgb = csfsurf(1, icsf) … … 5234 5241 IF ( isurfsrc == -1 ) THEN 5235 5242 ! 5236 !-- Diffuse rad from sky.5243 !-- Diffuse radiation from sky 5237 5244 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 5238 5245 ! 5239 !-- Absorbed diffuse LW from sky minus emitted to sky5246 !-- Absorbed diffuse LW radiation from sky minus emitted to sky 5240 5247 IF ( plant_lw_interact ) THEN 5241 5248 pcbinlw(ipcgb) = csf(1,icsf) & … … 5244 5251 ENDIF 5245 5252 ! 5246 !-- Direct rad5247 IF ( zenith(0)> 0 ) THEN5253 !-- Direct solar radiation 5254 IF ( cos_zenith > 0 ) THEN 5248 5255 !-- Estimate directed box absorption 5249 5256 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) … … 5312 5319 5313 5320 surfoutsl = albedo_surf * surfins 5314 !-- for non-transparent surfaces, longwave albedo is 1 - emissivity 5321 ! 5322 !-- for non-transparent surfaces, longwave albedo is 1 - emissivity 5315 5323 surfoutll = (1._wp - emiss_surf) * surfinl 5316 5324 … … 5318 5326 CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, & 5319 5327 surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 5320 IF ( ierr /= 0 ) THEN5328 IF ( ierr /= 0 ) THEN 5321 5329 WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, & 5322 5330 SIZE(surfouts), nsurfs, surfstart … … 5326 5334 CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, & 5327 5335 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) 5328 IF ( ierr /= 0 ) THEN5336 IF ( ierr /= 0 ) THEN 5329 5337 WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, & 5330 5338 SIZE(surfoutl), nsurfs, surfstart … … 5336 5344 surfoutl = surfoutll 5337 5345 #endif 5338 5339 !-- reset for next pass input5346 ! 5347 !-- Reset for the input from next reflective pass 5340 5348 surfins = 0._wp 5341 5349 surfinl = 0._wp 5342 5343 !-- reflected radiation5350 ! 5351 !-- Reflected radiation 5344 5352 DO isvf = 1, nsvfl 5345 5353 isurf = svfsurf(1, isvf) … … 5825 5833 !> This subroutine splits direct and diffusion dw radiation 5826 5834 !> It sould not be called in case the radiation model already does it 5827 !> It follows <CITATION>5835 !> It follows Boland, Ridley & Brown (2008) 5828 5836 !------------------------------------------------------------------------------! 5829 5837 SUBROUTINE calc_diffusion_radiation … … 5855 5863 !-- towards 0 while keeping full continuity. 5856 5864 !-- 5857 IF ( zenith(0)<= lowest_solarUp ) THEN5865 IF ( cos_zenith <= lowest_solarUp ) THEN 5858 5866 corrected_solarUp = lowest_solarUp 5859 5867 ELSE 5860 corrected_solarUp = zenith(0)5868 corrected_solarUp = cos_zenith 5861 5869 ENDIF 5862 5870 … … 8225 8233 !-- Update apparent solar position based on modified t_s_r_p 8226 8234 CALL calc_zenith 8227 IF ( zenith(0)> 0 ) THEN8235 IF ( cos_zenith > 0 ) THEN 8228 8236 !-- 8229 8237 !-- Identify solar direction vector (discretized number) 1) 8230 i = MODULO(NINT(ATAN2(sun_dir_lon (0), sun_dir_lat(0)) &8238 i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat) & 8231 8239 / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), & 8232 8240 raytrace_discrete_azims) 8233 j = FLOOR(ACOS( zenith(0)) / pi * raytrace_discrete_elevs)8241 j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs) 8234 8242 IF ( dsidir_rev(j, i) == -1 ) THEN 8235 8243 ndsidir = ndsidir + 1 -
palm/trunk/SOURCE/urban_surface_mod.f90
r3814 r3824 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Remove unused imports 31 ! 32 ! 33 ! 3814 2019-03-26 08:40:31Z pavelkrc 30 34 ! unused subroutine commented out 31 35 ! … … 438 442 439 443 USE basic_constants_and_equations_mod, & 440 ONLY: c_p, g, kappa, pi, r_d, rho_l, l_v 444 ONLY: c_p, g, kappa, pi, r_d, rho_l, l_v, sigma_sb 441 445 442 446 USE control_parameters, & … … 473 477 474 478 USE radiation_model_mod, & 475 ONLY: albedo_type, radiation_interaction, calc_zenith, zenith,&479 ONLY: albedo_type, radiation_interaction, & 476 480 radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, & 477 sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon, &478 481 force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u, & 479 482 iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id, &
Note: See TracChangeset
for help on using the changeset viewer.
