Changeset 3824 for palm/trunk/SOURCE/radiation_model_mod.f90
 Timestamp:
 Mar 27, 2019 3:56:16 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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 20190326 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 zcoordinate of solar unit vector 754 REAL(wp) :: sun_dir_lat !< ycoordinate of solar unit vector 755 REAL(wp) :: sun_dir_lon !< xcoordinate 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 nontransparent surfaces, longwave albedo is 1  emissivity 5321 ! 5322 ! for nontransparent 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
Note: See TracChangeset
for help on using the changeset viewer.