Changeset 3449 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- Oct 29, 2018 7:36:56 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (46 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r3435 r3449 28 28 ! ----------------- 29 29 ! $Id$ 30 ! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague) 31 ! - Interaction of plant canopy with LW radiation 32 ! - Transpiration from resolved plant canopy dependent on radiation 33 ! called from RTM 34 ! 35 ! 36 ! 3435 2018-10-26 18:25:44Z gronemeier 30 37 ! - workaround: return unit=illegal in check_data_output for certain variables 31 38 ! when check called from init_masks … … 499 506 500 507 USE basic_constants_and_equations_mod, & 501 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, 508 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, & 502 509 barometric_formula 503 510 … … 544 551 545 552 USE plant_canopy_model_mod, & 546 ONLY: lad_s, pc_heating_rate, pc_transpiration_rate 553 ONLY: lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate, & 554 plant_canopy_transpiration, pcm_calc_transpiration_rate 547 555 548 556 USE pegrid … … 853 861 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices 854 862 INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices 855 !< parameter but set in the code 863 REAL(wp), DIMENSION(0:nsurf_type) :: facearea !< area of single face in respective 864 !< direction (will be calc'd) 856 865 857 866 … … 887 896 !-- configuration parameters (they can be setup in PALM config) 888 897 LOGICAL :: raytrace_mpi_rma = .TRUE. !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing 889 LOGICAL :: read_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be read from a file at initialization 890 LOGICAL :: write_svf_on_init = .FALSE. !< flag parameter indicating wheather SVFs will be written out to a file 891 LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 898 LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 892 899 !< reflected radiation (as opposed to all mutually visible pairs) 900 LOGICAL :: plant_lw_interact = .TRUE. !< whether plant canopy interacts with LW radiation (in addition to SW) 893 901 INTEGER(iwp) :: mrt_nlevels = 0 !< number of vertical boxes above surface for which to calculate MRT 894 902 LOGICAL :: mrt_skip_roof = .TRUE. !< do not calculate MRT above roof surfaces 895 903 LOGICAL :: mrt_include_sw = .TRUE. !< should MRT calculation include SW radiation as well? 896 INTEGER(iwp) :: nrefsteps = 0!< number of reflection steps to perform904 INTEGER(iwp) :: nrefsteps = 3 !< number of reflection steps to perform 897 905 REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha) 898 906 INTEGER(iwp), PARAMETER :: rad_version_len = 10 !< length of identification string of rad version 899 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 1.1' !< identification of version of binary svf and restart files907 CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 3.0' !< identification of version of binary svf and restart files 900 908 INTEGER(iwp) :: raytrace_discrete_elevs = 40 !< number of discretization steps for elevation (nadir to zenith) 901 909 INTEGER(iwp) :: raytrace_discrete_azims = 80 !< number of discretization steps for azimuth (out of 360 degrees) … … 929 937 INTEGER(iwp) :: ity !< 930 938 INTEGER(iwp) :: itz !< 931 INTEGER(iwp) :: isurfs !< 932 REAL(wp) :: rsvf !< 939 INTEGER(iwp) :: isurfs !< Idx of source face / -1 for sky 940 REAL(wp) :: rcvf !< Canopy view factor for faces / 941 !< canopy sink factor for sky (-1) 933 942 END TYPE 934 943 … … 976 985 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection 977 986 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection 987 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlg !< global array of incoming lw radiation from plant canopy 978 988 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection 979 989 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection … … 2928 2938 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 2929 2939 2930 NAMELIST /radiation_par/ albedo, albedo_ type, albedo_lw_dir,&2931 albedo_ lw_dif, albedo_sw_dir, albedo_sw_dif,&2932 constant_albedo, dt_radiation, emissivity, &2933 lw_radiation, m rt_nlevels, mrt_skip_roof,&2934 m rt_include_sw, net_radiation,&2935 radiation_scheme, skip_time_do_radiation,&2936 sw_radiation, unscheduled_radiation_calls,&2937 max_raytracing_dist, min_irrf_value,&2938 nrefsteps, raytrace_mpi_rma,&2939 surface_reflections, svfnorm_report_thresh,&2940 radiation_interactions_on,&2941 rad_angular_discretization,&2942 raytrace_discrete_azims, &2943 raytrace_discrete_elevs 2940 NAMELIST /radiation_par/ albedo, albedo_lw_dif, albedo_lw_dir, & 2941 albedo_sw_dif, albedo_sw_dir, albedo_type, & 2942 constant_albedo, dt_radiation, emissivity, & 2943 lw_radiation, max_raytracing_dist, & 2944 min_irrf_value, mrt_include_sw, mrt_nlevels, & 2945 mrt_skip_roof, net_radiation, nrefsteps, & 2946 plant_lw_interact, rad_angular_discretization,& 2947 radiation_interactions_on, radiation_scheme, & 2948 raytrace_discrete_azims, & 2949 raytrace_discrete_elevs, raytrace_mpi_rma, & 2950 skip_time_do_radiation, surface_reflections, & 2951 svfnorm_report_thresh, sw_radiation, & 2952 unscheduled_radiation_calls 2953 2944 2954 2945 NAMELIST /radiation_parameters/ albedo, albedo_type, albedo_lw_dir, & 2946 albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, & 2947 constant_albedo, dt_radiation, emissivity, & 2948 lw_radiation, mrt_nlevels, mrt_skip_roof, & 2949 mrt_include_sw, net_radiation, & 2950 radiation_scheme, skip_time_do_radiation, & 2951 sw_radiation, unscheduled_radiation_calls, & 2952 max_raytracing_dist, min_irrf_value, & 2953 nrefsteps, raytrace_mpi_rma, & 2954 surface_reflections, svfnorm_report_thresh, & 2955 radiation_interactions_on, & 2956 rad_angular_discretization, & 2957 raytrace_discrete_azims, & 2958 raytrace_discrete_elevs 2955 NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir, & 2956 albedo_sw_dif, albedo_sw_dir, albedo_type, & 2957 constant_albedo, dt_radiation, emissivity, & 2958 lw_radiation, max_raytracing_dist, & 2959 min_irrf_value, mrt_include_sw, mrt_nlevels, & 2960 mrt_skip_roof, net_radiation, nrefsteps, & 2961 plant_lw_interact, rad_angular_discretization,& 2962 radiation_interactions_on, radiation_scheme, & 2963 raytrace_discrete_azims, & 2964 raytrace_discrete_elevs, raytrace_mpi_rma, & 2965 skip_time_do_radiation, surface_reflections, & 2966 svfnorm_report_thresh, sw_radiation, & 2967 unscheduled_radiation_calls 2959 2968 2960 2969 line = ' ' … … 4648 4657 REAL(wp), DIMENSION(0:nsurf_type) :: costheta !< direct irradiance factor of solar angle 4649 4658 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temperature tendency 4650 REAL(wp), DIMENSION(nzub:nzut) :: pctf_prep !< precalculated factor for canopy transpiration tendency4651 4659 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4652 4660 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 4653 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 4661 REAL(wp) :: asrc !< area of source face 4662 REAL(wp) :: pcrad !< irradiance from plant canopy 4654 4663 REAL(wp) :: pabsswl = 0.0_wp !< total absorbed SW radiation energy in local processor (W) 4655 4664 REAL(wp) :: pabssw = 0.0_wp !< total absorbed SW radiation energy in all processors (W) … … 4672 4681 pchf_prep(:) = r_d * exner(nzub:nzut) & 4673 4682 / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T) 4674 pctf_prep(:) = r_d * exner(nzub:nzut) &4675 / (l_v * hyp(nzub:nzut) * dx*dy*dz(1))4676 4683 ENDIF 4677 4684 #endif … … 4729 4736 mrtinlw(:) = 0._wp 4730 4737 ENDIF 4738 surfinlg(:) = 0._wp !global 4731 4739 4732 4740 … … 4824 4832 ! 4825 4833 !-- For surface-to-surface factors we calculate thermal radiation in 1st pass 4826 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4834 IF ( plant_lw_interact ) THEN 4835 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc) 4836 ELSE 4837 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4838 ENDIF 4827 4839 ENDDO 4828 4840 ENDIF … … 4833 4845 i = surfl(ix, isurf) 4834 4846 surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf) 4835 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 4847 IF ( plant_lw_interact ) THEN 4848 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf) 4849 ELSE 4850 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 4851 ENDIF 4836 4852 ENDDO 4837 4853 ! … … 4880 4896 pcbinswdir(:) = 0._wp 4881 4897 pcbinswdif(:) = 0._wp 4882 pcbinlw(:) = 0._wp !< will stay always 0 since we don't absorb lw anymore4883 ! 4884 !-- 4898 pcbinlw(:) = 0._wp 4899 ! 4900 !-- pcsf first pass 4885 4901 DO icsf = 1, ncsfl 4886 4902 ipcgb = csfsurf(1, icsf) … … 4891 4907 4892 4908 IF ( isurfsrc == -1 ) THEN 4893 !-- Diffuse rad from sky. 4894 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 4895 4896 !--Direct rad 4897 IF ( zenith(0) > 0 ) THEN 4898 !--Estimate directed box absorption 4899 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) 4900 4901 !--isd has already been established, see 1) 4902 pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area & 4903 * pc_abs_frac * dsitransc(ipcgb, isd) 4904 ENDIF 4909 ! 4910 !-- Diffuse rad from sky. 4911 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 4912 ! 4913 !-- Absorbed diffuse LW from sky minus emitted to sky 4914 IF ( plant_lw_interact ) THEN 4915 pcbinlw(ipcgb) = csf(1,icsf) & 4916 * (rad_lw_in_diff(j, i) & 4917 - sigma_sb * (pt(k,j,i)*exner(k))**4) 4918 ENDIF 4919 ! 4920 !-- Direct rad 4921 IF ( zenith(0) > 0 ) THEN 4922 !-- Estimate directed box absorption 4923 pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i)) 4924 ! 4925 !-- isd has already been established, see 1) 4926 pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area & 4927 * pc_abs_frac * dsitransc(ipcgb, isd) 4928 ENDIF 4929 ELSE 4930 IF ( plant_lw_interact ) THEN 4931 ! 4932 !-- Thermal emission from plan canopy towards respective face 4933 pcrad = sigma_sb * (pt(k,j,i) * exner(k))**4 * csf(1,icsf) 4934 surfinlg(isurfsrc) = surfinlg(isurfsrc) + pcrad 4935 ! 4936 !-- Remove the flux above + absorb LW from first pass from surfaces 4937 asrc = facearea(surf(id, isurfsrc)) 4938 pcbinlw(ipcgb) = pcbinlw(ipcgb) & 4939 + (csf(1,icsf) * surfoutl(isurfsrc) & ! Absorb from first pass surf emit 4940 - pcrad) & ! Remove emitted heatflux 4941 * asrc 4942 ENDIF 4905 4943 ENDIF 4906 4944 ENDDO … … 4908 4946 pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:) 4909 4947 ENDIF 4948 4949 IF ( plant_lw_interact ) THEN 4950 ! 4951 !-- Exchange incoming lw radiation from plant canopy 4952 #if defined( __parallel ) 4953 CALL MPI_Allreduce(MPI_IN_PLACE, surfinlg, nsurf, MPI_REAL, MPI_SUM, comm2d, ierr) 4954 IF ( ierr /= 0 ) THEN 4955 WRITE (9,*) 'Error MPI_Allreduce:', ierr 4956 FLUSH(9) 4957 ENDIF 4958 surfinl(:) = surfinl(:) + surfinlg(surfstart(myid)+1:surfstart(myid+1)) 4959 #else 4960 surfinl(:) = surfinl(:) + surfinlg(:) 4961 #endif 4962 ENDIF 4963 4910 4964 surfins = surfinswdir + surfinswdif 4911 4965 surfinl = surfinl + surfinlwdif … … 4966 5020 isurfsrc = svfsurf(2, isvf) 4967 5021 surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc) 4968 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 5022 IF ( plant_lw_interact ) THEN 5023 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc) 5024 ELSE 5025 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 5026 ENDIF 4969 5027 ENDDO 4970 4971 !-- radiation absorbed by plant canopy 4972 DO icsf = 1, ncsfl 5028 ! 5029 !-- NOTE: PC absorbtion and MRT from reflected can both be done at once 5030 !-- after all reflections if we do one more MPI_ALLGATHERV on surfout. 5031 !-- Advantage: less local computation. Disadvantage: one more collective 5032 !-- MPI call. 5033 ! 5034 !-- Radiation absorbed by plant canopy 5035 DO icsf = 1, ncsfl 4973 5036 ipcgb = csfsurf(1, icsf) 4974 5037 isurfsrc = csfsurf(2, icsf) 4975 5038 IF ( isurfsrc == -1 ) CYCLE ! sky->face only in 1st pass, not here 4976 4977 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) 5039 ! 5040 !-- Calculate source surface area. If the `surf' array is removed 5041 !-- before timestepping starts (future version), then asrc must be 5042 !-- stored within `csf' 5043 asrc = facearea(surf(id, isurfsrc)) 5044 pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) * asrc 5045 IF ( plant_lw_interact ) THEN 5046 pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc) * asrc 5047 ENDIF 4978 5048 ENDDO 4979 5049 ! … … 4996 5066 IF ( npcbl > 0 ) THEN 4997 5067 pc_heating_rate(:,:,:) = 0.0_wp 4998 pc_transpiration_rate(:,:,:) = 0.0_wp4999 5068 DO ipcgb = 1, npcbl 5000 5001 5069 j = pcbl(iy, ipcgb) 5002 5070 i = pcbl(ix, ipcgb) 5003 5071 k = pcbl(iz, ipcgb) 5004 5072 ! 5005 !-- 5073 !-- Following expression equals former kk = k - nzb_s_inner(j,i) 5006 5074 kk = k - get_topography_top_index_ji( j, i, 's' ) !- lad arrays are defined flat 5007 5075 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 5008 5076 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt 5009 5010 ! pc_transpiration_rate(kk,j,i) = 0.75_wp* (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &5011 ! * pctf_prep(k) * pt(k, j, i) !-- = dq/dt5012 5013 5077 ENDDO 5078 5079 IF ( humidity .AND. plant_canopy_transpiration ) THEN 5080 !-- Calculation of plant canopy transpiration rate and correspondidng latent heat rate 5081 pc_transpiration_rate(:,:,:) = 0.0_wp 5082 pc_latent_rate(:,:,:) = 0.0_wp 5083 DO ipcgb = 1, npcbl 5084 i = pcbl(ix, ipcgb) 5085 j = pcbl(iy, ipcgb) 5086 k = pcbl(iz, ipcgb) 5087 kk = k - get_topography_top_index_ji( j, i, 's' ) !- lad arrays are defined flat 5088 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), & 5089 pc_transpiration_rate(kk,j,i), pc_latent_rate(kk,j,i) ) 5090 ENDDO 5091 ENDIF 5014 5092 ENDIF 5015 5093 ! … … 5161 5239 !-- domain when using average_radiation in the respective radiation model 5162 5240 5163 !-- Precalculate face areas for all face directions using normal vector5164 DO d = 0, nsurf_type5165 facearea(d) = 1._wp5166 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx5167 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy5168 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)5169 ENDDO5170 5241 !-- calculate horizontal area 5171 5242 ! !!! ATTENTION!!! uniform grid is assumed here … … 5429 5500 IMPLICIT NONE 5430 5501 5431 INTEGER(iwp) :: i, j, k, l, m 5502 INTEGER(iwp) :: i, j, k, l, m, d 5432 5503 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 5433 5504 INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt … … 5439 5510 #endif 5440 5511 5512 ! 5513 !-- precalculate face areas for different face directions using normal vector 5514 DO d = 0, nsurf_type 5515 facearea(d) = 1._wp 5516 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 5517 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 5518 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1) 5519 ENDDO 5441 5520 ! 5442 5521 !-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be … … 5904 5983 ALLOCATE( surfouts(nsurf) ) 5905 5984 ALLOCATE( surfoutl(nsurf) ) 5985 ALLOCATE( surfinlg(nsurf) ) 5906 5986 ALLOCATE( skyvf(nsurfl) ) 5907 5987 ALLOCATE( skyvft(nsurfl) ) … … 5940 6020 REAL(wp) :: az1, az2 !< relative azimuth of section borders 5941 6021 REAL(wp) :: azmid !< ray (center) azimuth 6022 REAL(wp) :: yxlen !< |yxdir| 5942 6023 REAL(wp), DIMENSION(2) :: yxdir !< y,x *unit* vector of ray direction (in grid units) 5943 6024 REAL(wp), DIMENSION(:), ALLOCATABLE :: zdirs !< directions in z (tangent of elevation) 6025 REAL(wp), DIMENSION(:), ALLOCATABLE :: zcent !< zenith angle centers 5944 6026 REAL(wp), DIMENSION(:), ALLOCATABLE :: zbdry !< zenith angle boundaries 5945 6027 REAL(wp), DIMENSION(:), ALLOCATABLE :: vffrac !< view factor fractions for individual rays … … 5949 6031 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: itarget !< face indices of detected obstacles 5950 6032 INTEGER(iwp) :: itarg0, itarg1 5951 #if defined( __parallel ) 5952 #endif 5953 5954 5955 5956 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 6033 5957 6034 INTEGER(iwp) :: udim 5958 6035 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l … … 5966 6043 LOGICAL :: visible 5967 6044 REAL(wp), DIMENSION(3) :: sa, ta !< real coordinates z,y,x of source and target 6045 REAL(wp) :: difvf !< differential view factor 5968 6046 REAL(wp) :: transparency, rirrf, sqdist, svfsum 5969 6047 INTEGER(iwp) :: isurflt, isurfs, isurflt_prev … … 5983 6061 CALL location_message( ' calculation of SVF and CSF', .TRUE. ) 5984 6062 CALL radiation_write_debug_log('Start calculation of SVF and CSF') 5985 5986 !-- precalculate face areas for different face directions using normal vector5987 DO d = 0, nsurf_type5988 facearea(d) = 1._wp5989 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx5990 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy5991 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)5992 ENDDO5993 6063 5994 6064 !-- initialize variables and temporary arrays for calculation of svf and csf … … 6219 6289 END SELECT 6220 6290 6221 ALLOCATE ( zdirs(1:nzn), z bdry(0:nzn), vffrac(1:nzn*naz), &6291 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), & 6222 6292 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6223 6293 !in case of rad_angular_discretization … … 6225 6295 itarg0 = 1 6226 6296 itarg1 = nzn 6227 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6297 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6228 6298 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 6229 6299 IF ( td == iup_u .OR. td == iup_l ) THEN … … 6251 6321 !-- sum of whole vffrac equals 1, verified 6252 6322 ENDIF 6253 yxdir = (/ COS(azmid), SIN(azmid) /) 6323 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6324 yxlen = SQRT(SUM(yxdir(:)**2)) 6325 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6326 yxdir(:) = yxdir(:) / yxlen 6327 6254 6328 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6255 6329 surfstart(myid) + isurflt, facearea(td), & … … 6257 6331 .FALSE., lowest_free_ray, & 6258 6332 ztransp(itarg0:itarg1), & 6259 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6260 !FIXME itarget available only in 6261 !case of rad_angular_discretization 6333 itarget(itarg0:itarg1)) 6334 6262 6335 skyvf(isurflt) = skyvf(isurflt) + & 6263 6336 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1)) … … 6337 6410 ENDIF ! rad_angular_discretization 6338 6411 6339 DEALLOCATE ( zdirs, z bdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only6412 DEALLOCATE ( zdirs, zcent, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only 6340 6413 !in case of rad_angular_discretization 6341 6414 ! … … 6367 6440 ENDIF 6368 6441 6442 difvf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 6443 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 6444 / (pi * sqdist) ! square of distance between centers 6445 ! 6369 6446 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 6370 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 6371 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 6372 / (pi * sqdist) & ! square of distance between centers 6373 * facearea(sd) 6447 rirrf = difvf * facearea(sd) 6374 6448 6375 6449 !-- reject raytracing for potentially too small view factor values … … 6380 6454 6381 6455 !-- raytrace + process plant canopy sinks within 6382 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &6456 CALL raytrace(sa, ta, isurfs, difvf, facearea(td), .TRUE., & 6383 6457 visible, transparency) 6384 6458 … … 6428 6502 nzn = raytrace_discrete_elevs / 2 6429 6503 zns = pi / 2._wp / REAL(nzn, wp) 6430 ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn), &6504 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), vffrac(1:nzn), ztransp(1:nzn), & 6431 6505 itarget(1:nzn) ) 6432 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6506 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6433 6507 vffrac(:) = 0._wp 6434 6508 … … 6440 6514 DO iaz = 1, naz 6441 6515 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 6442 yxdir = (/ COS(azmid), SIN(azmid) /) 6516 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6517 yxlen = SQRT(SUM(yxdir(:)**2)) 6518 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6519 yxdir(:) = yxdir(:) / yxlen 6443 6520 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6444 6521 -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., & 6445 lowest_free_ray, ztransp, itarget) !FIXME unit vect in grid units + zdirs6522 lowest_free_ray, ztransp, itarget) 6446 6523 6447 6524 !-- Save direct solar transparency … … 6456 6533 ENDDO 6457 6534 ENDDO 6458 DEALLOCATE ( zdirs, vffrac, ztransp, itarget )6535 DEALLOCATE ( zdirs, zcent, vffrac, ztransp, itarget ) 6459 6536 !-- 6460 6537 !-- Raytrace to MRT boxes … … 6469 6546 nzn = raytrace_discrete_elevs 6470 6547 zns = pi / REAL(nzn, wp) 6471 ALLOCATE ( zdirs(1:nzn), z bdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &6548 ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), & 6472 6549 ztransp(1:nzn*naz), itarget(1:nzn*naz) ) !FIXME allocate itarget only 6473 6550 !in case of rad_angular_discretization 6474 6551 6475 z dirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)6552 zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/) 6476 6553 zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/) 6477 6554 vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp) … … 6493 6570 DO iaz = 1, naz 6494 6571 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs 6495 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), nzn, zdirs, & 6572 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /) 6573 yxlen = SQRT(SUM(yxdir(:)**2)) 6574 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:))) 6575 yxdir(:) = yxdir(:) / yxlen 6576 6577 CALL raytrace_2d(ta, yxdir, nzn, zdirs, & 6496 6578 -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., & 6497 6579 .FALSE., .TRUE., lowest_free_ray, & 6498 6580 ztransp(itarg0:itarg1), & 6499 itarget(itarg0:itarg1)) !FIXME unit vect in grid units + zdirs 6500 !FIXME itarget available only in 6501 !case of rad_angular_discretization 6581 itarget(itarg0:itarg1)) 6502 6582 6503 6583 !-- Sky view factors for MRT … … 6574 6654 6575 6655 ENDDO ! imrt 6576 DEALLOCATE ( zdirs, z bdry, vffrac, vffrac0, ztransp, itarget )6656 DEALLOCATE ( zdirs, zcent, zbdry, vffrac, vffrac0, ztransp, itarget ) 6577 6657 ! 6578 6658 !-- Move MRT factors to final arrays … … 6778 6858 j = 0 6779 6859 ENDIF 6780 !-- fill out real values of rsvf, rtransp 6781 csflt(1,kcsf) = acsf(kcsf)%rsvf 6860 csflt(1,kcsf) = acsf(kcsf)%rcvf 6782 6861 !-- fill out integer values of itz,ity,itx,isurfs 6783 6862 kcsflt(1,kcsf) = acsf(kcsf)%itz … … 6946 7025 !> doesn't need to be the same in all three dimensions). 6947 7026 !------------------------------------------------------------------------------! 6948 SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency)7027 SUBROUTINE raytrace(src, targ, isrc, difvf, atarg, create_csf, visible, transparency) 6949 7028 IMPLICIT NONE 6950 7029 6951 7030 REAL(wp), DIMENSION(3), INTENT(in) :: src, targ !< real coordinates z,y,x 6952 7031 INTEGER(iwp), INTENT(in) :: isrc !< index of source face for csf 6953 REAL(wp), INTENT(in) :: rirrf !< irradiancefactor for csf7032 REAL(wp), INTENT(in) :: difvf !< differential view factor for csf 6954 7033 REAL(wp), INTENT(in) :: atarg !< target surface area for csf 6955 7034 LOGICAL, INTENT(in) :: create_csf !< whether to generate new CSFs during raytracing … … 7113 7192 acsf(ncsfl)%itz = boxes(1,i) 7114 7193 acsf(ncsfl)%isurfs = isrc 7115 acsf(ncsfl)%r svf = cursink*transparency*rirrf*atarg7194 acsf(ncsfl)%rcvf = cursink*transparency*difvf*atarg 7116 7195 ENDIF !< create_csf 7117 7196 … … 7514 7593 acsf(ncsfl)%itz = iz 7515 7594 acsf(ncsfl)%isurfs = iorig 7516 acsf(ncsfl)%r svf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)7595 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k) 7517 7596 ENDIF 7518 7597 … … 7571 7650 acsf(ncsfl)%ity = rt2_track(1,i) 7572 7651 acsf(ncsfl)%itz = iz 7573 acsf(ncsfl)%isurfs = itarget(k) !if above horizon, then itarget(k)==-1, which7574 !is also a special ID indicating sky7575 acsf(ncsfl)%r svf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)7652 IF ( itarget(k) /= -1 ) ERROR STOP !FIXME remove after test 7653 acsf(ncsfl)%isurfs = -1 7654 acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k) 7576 7655 ENDIF ! create_csf 7577 7656 … … 7600 7679 INTEGER(iwp), TARGET, INTENT(out) :: isurfl 7601 7680 INTEGER(iwp), INTENT(out) :: iproc 7681 #if defined( __parallel ) 7682 #else 7683 INTEGER(iwp) :: target_displ !< index of the grid in the local gridsurf array 7684 #endif 7602 7685 INTEGER(iwp) :: px, py !< number of processors in x and y direction 7603 7686 !< before the processor in the question … … 8138 8221 .AND. acsfnew(iwrite)%isurfs == acsf(iread)%isurfs ) THEN 8139 8222 8140 acsfnew(iwrite)%r svf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf8223 acsfnew(iwrite)%rcvf = acsfnew(iwrite)%rcvf + acsf(iread)%rcvf 8141 8224 !-- advance reading index, keep writing index 8142 8225 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
