Changeset 4529 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- May 12, 2020 9:14:57 AM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r4517 r4529 28 28 ! ----------------- 29 29 ! $Id$ 30 ! - added the following new features to the coupling of RTM-radiation model: 31 ! 1) considering the vegetation interaction with LW in the coupling 32 ! 2) considering PC emissivity in calculating the effective emissivity 33 ! 3) new algorithm for claculating the coupling parameters so that each term 34 ! is calculated within its original line and not at the end. 35 ! - minor formatting and comments changes 36 ! 37 ! 4517 2020-05-03 14:29:30Z raasch 30 38 ! added restart with MPI-IO for reading local arrays 31 39 ! … … 5805 5813 REAL(wp) :: asrc !< area of source face 5806 5814 REAL(wp) :: pcrad !< irradiance from plant canopy 5807 REAL(wp) :: pabsswl = 0.0_wp !< total absorbed SW radiation energy in local processor (W) 5808 REAL(wp) :: pabssw = 0.0_wp !< total absorbed SW radiation energy in all processors (W) 5809 REAL(wp) :: pabslwl = 0.0_wp !< total absorbed LW radiation energy in local processor (W) 5810 REAL(wp) :: pabslw = 0.0_wp !< total absorbed LW radiation energy in all processors (W) 5811 REAL(wp) :: pemitlwl = 0.0_wp !< total emitted LW radiation energy in all processors (W) 5812 REAL(wp) :: pemitlw = 0.0_wp !< total emitted LW radiation energy in all processors (W) 5813 REAL(wp) :: pinswl = 0.0_wp !< total received SW radiation energy in local processor (W) 5814 REAL(wp) :: pinsw = 0.0_wp !< total received SW radiation energy in all processor (W) 5815 REAL(wp) :: pinlwl = 0.0_wp !< total received LW radiation energy in local processor (W) 5816 REAL(wp) :: pinlw = 0.0_wp !< total received LW radiation energy in all processor (W) 5817 REAL(wp) :: emiss_sum_surfl !< sum of emissisivity of surfaces in local processor 5818 REAL(wp) :: emiss_sum_surf !< sum of emissisivity of surfaces in all processor 5819 REAL(wp) :: area_surfl !< total area of surfaces in local processor 5820 REAL(wp) :: area_surf !< total area of surfaces in all processor 5821 REAL(wp) :: area_hor !< total horizontal area of domain in all processor 5815 !- variables for coupling the radiation modle (e.g. RRTMG) and RTM 5816 REAL(wp) :: pabsswl !< total absorbed SW radiation energy in local processor (W) 5817 REAL(wp) :: pabssw !< total absorbed SW radiation energy in all processors (W) 5818 REAL(wp) :: pabslwl !< total absorbed LW radiation energy in local processor (W) 5819 REAL(wp) :: pabslw !< total absorbed LW radiation energy in all processors (W) 5820 REAL(wp) :: pemitlwl !< total emitted LW radiation energy in all processors (W) 5821 REAL(wp) :: pemitlw !< total emitted LW radiation energy in all processors (W) 5822 REAL(wp) :: pinswl !< total received SW radiation energy in local processor (W) 5823 REAL(wp) :: pinsw !< total received SW radiation energy in all processor (W) 5824 REAL(wp) :: pinlwl !< total received LW radiation energy in local processor (W) 5825 REAL(wp) :: pinlw !< total received LW radiation energy in all processor (W) 5826 REAL(wp) :: area_norm !< reference horizontal area of domain in all processor 5827 REAL(wp) :: pabs_surf_lwdifl !< total absorbed LW radiation in surfaces from sky in local processor (W) 5828 REAL(wp) :: pabs_surf_lwdif !< total absorbed LW radiation in surfaces from sky in all processors (W) 5829 REAL(wp) :: pabs_pc_lwdifl !< total absorbed LW radiation in plant canopy from sky in local processor (W) 5830 REAL(wp) :: pabs_pc_lwdif !< total absorbed LW radiation in plant canopy from sky in all processors (W) 5831 5822 5832 REAL(wp) :: sun_direct_factor !< factor for direct normal radiation from direct horizontal 5823 5833 #if defined( __parallel ) … … 5888 5898 ENDIF 5889 5899 5890 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5891 !-- First pass: direct + diffuse irradiance + thermal 5892 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5893 surfinswdir = 0._wp !nsurfl 5894 surfins = 0._wp !nsurfl 5895 surfinl = 0._wp !nsurfl 5896 surfoutsl(:) = 0.0_wp !start-end 5897 surfoutll(:) = 0.0_wp !start-end 5900 5901 !-- First pass of radiation interaction: 5902 !-- 1) direct and diffuse irradiance 5903 !-- 2) thermal emissions 5904 5905 !- Initialize relavant surface flux arrays and radiation energy sum 5906 !- surface flux 5907 surfinswdir = 0.0_wp 5908 surfins = 0.0_wp 5909 surfinl = 0.0_wp 5910 surfoutsl(:) = 0.0_wp 5911 surfoutll(:) = 0.0_wp 5898 5912 IF ( nmrtbl > 0 ) THEN 5899 mrtinsw(:) = 0. _wp5900 mrtinlw(:) = 0. _wp5913 mrtinsw(:) = 0.0_wp 5914 mrtinlw(:) = 0.0_wp 5901 5915 ENDIF 5902 surfinlg(:) = 0._wp !global 5903 5916 surfinlg(:) = 0.0_wp 5917 !- radiation energy sum 5918 pinlwl = 0.0_wp 5919 pinswl = 0.0_wp 5920 pemitlwl = 0.0_wp 5921 pabsswl = 0.0_wp 5922 pabslwl = 0.0_wp 5923 pabs_surf_lwdifl = 0.0_wp 5924 pabs_pc_lwdifl = 0.0_wp 5904 5925 5905 5926 !-- Set up thermal radiation from surfaces … … 6019 6040 j = surfl(iy, isurf) 6020 6041 i = surfl(ix, isurf) 6042 d = surfl(id, isurf) 6021 6043 surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf) 6044 !- update received SW energy for RTM coupling 6045 pinswl = pinswl + surfinswdif(isurf) * facearea(d) 6022 6046 IF ( plant_lw_interact ) THEN 6023 6047 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf) 6048 !- update received LW energy for RTM coupling 6049 pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d) 6024 6050 ELSE 6025 6051 surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf) 6052 !- update received LW energy for RTM coupling 6053 pinlwl = pinlwl + surfinlwdif(isurf) * facearea(d) 6026 6054 ENDIF 6027 6055 ENDDO … … 6053 6081 j = surfl(iy, isurf) 6054 6082 i = surfl(ix, isurf) 6083 d = surfl(id, isurf) 6055 6084 surfinswdir(isurf) = rad_sw_in_dir(j,i) * & 6056 6085 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) * sun_direct_factor 6086 !- update received SW energy for RTM coupling 6087 pinswl = pinswl + surfinswdir(isurf) * facearea(d) 6057 6088 ENDDO 6058 6089 ! … … 6092 6123 !-- Diffuse radiation from sky 6093 6124 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i) 6125 !-- add to the sum of SW radiation energy 6126 pinswl = pinswl + pcbinswdif(ipcgb) 6094 6127 ! 6095 6128 !-- Absorbed diffuse LW radiation from sky minus emitted to sky … … 6098 6131 * (rad_lw_in_diff(j, i) & 6099 6132 - sigma_sb * (pt(k,j,i)*exner(k))**4) 6133 pinlwl = pinlwl + csf(1,icsf) * rad_lw_in_diff(j,i) 6134 pabslwl = pabslwl + csf(1,icsf) * rad_lw_in_diff(j,i) 6135 pemitlwl = pemitlwl + & 6136 csf(1,icsf) * sigma_sb * (pt(k,j,i)*exner(k))**4 6137 pabs_pc_lwdifl = pabs_pc_lwdifl + & 6138 csf(1,icsf) * rad_lw_in_diff(j,i) 6100 6139 ENDIF 6101 6140 ! … … 6108 6147 pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area & 6109 6148 * pc_abs_frac * dsitransc(ipcgb, isd) 6149 !-- add to the sum of SW radiation energy 6150 pinswl = pinswl + pcbinswdir(ipcgb) 6110 6151 ENDIF 6111 6152 ELSE … … 6122 6163 - pcrad) & ! Remove emitted heatflux 6123 6164 * asrc 6165 pabslwl = pabslwl + csf(1,icsf) * surfoutl(isurfsrc) * asrc 6166 pemitlwl = pemitlwl + pcrad * asrc 6124 6167 ENDIF 6125 6168 ENDIF … … 6178 6221 ENDIF 6179 6222 6180 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!6181 !-- Next passes -reflections6182 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 6223 !-- Next passes of radiation interactions: 6224 !-- radiation reflections 6225 6183 6226 DO refstep = 1, nrefsteps 6184 6227 … … 6273 6316 6274 6317 ENDDO ! refstep 6275 6276 !-- push heat flux absorbed by plant canopy to respective 3D arrays 6318 ! 6319 !-- push heat flux absorbed by plant canopy to respective 3D arrays and 6320 !-- add absorbed SW radiation energy for RTM coupling variables 6277 6321 IF ( npcbl > 0 ) THEN 6278 6322 pcm_heating_rate(:,:,:) = 0.0_wp … … 6286 6330 pcm_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 6287 6331 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt 6332 !-- add the absorbed SW radiation energy by plant canopy 6333 pabsswl = pabsswl + pcbinsw(ipcgb) 6288 6334 ENDDO 6289 6335 … … 6297 6343 k = pcbl(iz, ipcgb) 6298 6344 kk = k - topo_top_ind(j,i,0) !- lad arrays are defined flat 6299 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), & 6300 pcm_transpiration_rate(kk,j,i), pcm_latent_rate(kk,j,i) ) 6345 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), & 6346 pcbinlw(ipcgb), & 6347 pcm_transpiration_rate(kk,j,i), & 6348 pcm_latent_rate(kk,j,i) ) 6301 6349 ENDDO 6302 6350 ENDIF … … 6312 6360 ENDIF 6313 6361 ! 6314 !-- Transfer radiation arrays required for energy balance to the respective data types 6362 !-- Transfer radiation arrays required for energy balance to the respective data types 6363 ! and claculate relevant radiation model-RTM coupling terms 6364 6315 6365 DO i = 1, nsurfl 6316 6366 m = surfl(im,i) 6367 d = surfl(id, i) 6317 6368 ! 6318 6369 !-- (1) Urban surfaces … … 6493 6544 surf_lsm_v(3)%rad_lw_res(m) = surfinl(i) 6494 6545 ENDIF 6495 6546 ! 6547 !-- RTM coupling terms 6548 !-- sum of absorbed SW & LW radiation energy 6549 pabsswl = pabsswl + & 6550 (1.0_wp - albedo_surf(i)) * surfinsw(i) * facearea(d) 6551 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d) 6552 !-- sum of emitted LW radiation energy 6553 pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d) 6554 !-- emiss1 6555 pabs_surf_lwdifl = pabs_surf_lwdifl + & 6556 emiss_surf(i) * facearea(d) * surfinlwdif(i) 6496 6557 ENDDO 6497 6558 … … 6527 6588 ENDDO 6528 6589 ! 6529 !-- Calculate the average temperature, albedo, and emissivity for urban/land6530 !-- domain when using average_radiation in the respective radiation model6531 6532 !-- calculate horizontal area6533 ! !!! ATTENTION!!! uniform grid is assumed here6534 area_hor = (nx+1) * (ny+1) * dx * dy6535 !6536 !-- absorbed/received SW & LW and emitted LW energy of all physical6537 !-- surfaces (land and urban) in local processor6538 pinswl = 0._wp6539 pinlwl = 0._wp6540 pabsswl = 0._wp6541 pabslwl = 0._wp6542 pemitlwl = 0._wp6543 emiss_sum_surfl = 0._wp6544 area_surfl = 0._wp6545 DO i = 1, nsurfl6546 d = surfl(id, i)6547 !-- received SW & LW6548 pinswl = pinswl + (surfinswdir(i) + surfinswdif(i)) * facearea(d)6549 pinlwl = pinlwl + surfinlwdif(i) * facearea(d)6550 !-- absorbed SW & LW6551 pabsswl = pabsswl + (1._wp - albedo_surf(i)) * &6552 surfinsw(i) * facearea(d)6553 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)6554 !-- emitted LW6555 pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)6556 !-- emissivity and area sum6557 emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)6558 area_surfl = area_surfl + facearea(d)6559 END DO6560 !6561 !-- add the absorbed SW energy by plant canopy6562 IF ( npcbl > 0 ) THEN6563 pabsswl = pabsswl + SUM(pcbinsw)6564 pabslwl = pabslwl + SUM(pcbinlw)6565 pinswl = pinswl + SUM(pcbinswdir) + SUM(pcbinswdif)6566 ENDIF6567 !6568 6590 !-- gather all rad flux energy in all processors. In order to reduce 6569 6591 !-- the number of MPI calls (to reduce latencies), combine the required … … 6576 6598 combine_allreduce_l(4) = pabslwl 6577 6599 combine_allreduce_l(5) = pemitlwl 6578 combine_allreduce_l(6) = emiss_sum_surfl6579 combine_allreduce_l(7) = area_surfl6600 combine_allreduce_l(6) = pabs_surf_lwdifl 6601 combine_allreduce_l(7) = pabs_pc_lwdifl 6580 6602 6581 6603 CALL MPI_ALLREDUCE( combine_allreduce_l, & … … 6587 6609 ierr ) 6588 6610 6589 pinsw = combine_allreduce(1)6590 pinlw = combine_allreduce(2)6591 pabssw = combine_allreduce(3)6592 pabslw = combine_allreduce(4)6593 pemitlw = combine_allreduce(5)6594 emiss_sum_surf = combine_allreduce(6)6595 area_surf= combine_allreduce(7)6611 pinsw = combine_allreduce(1) 6612 pinlw = combine_allreduce(2) 6613 pabssw = combine_allreduce(3) 6614 pabslw = combine_allreduce(4) 6615 pemitlw = combine_allreduce(5) 6616 pabs_surf_lwdif = combine_allreduce(6) 6617 pabs_pc_lwdif = combine_allreduce(7) 6596 6618 #else 6597 pinsw = pinswl6598 pinlw = pinlwl6599 pabssw = pabsswl6600 pabslw = pabslwl6601 pemitlw = pemitlwl6602 emiss_sum_surf = emiss_sum_surfl6603 area_surf = area_surfl6619 pinsw = pinswl 6620 pinlw = pinlwl 6621 pabssw = pabsswl 6622 pabslw = pabslwl 6623 pemitlw = pemitlwl 6624 pabs_surf_lwdif = pabs_surf_lwdifl 6625 pabs_pc_lwdif = pabs_pc_lwdif 6604 6626 #endif 6605 6606 !-- (1) albedo 6627 ! 6628 !-- Calculate the effective radiation surface parameters based on 6629 !-- the parameterizations in Krc et al. 2020 6630 6631 !- (1) albedo Eq. * in Krc et al. 2020 6607 6632 IF ( pinsw /= 0.0_wp ) albedo_urb = ( pinsw - pabssw ) / pinsw 6608 !-- (2) average emmsivity 6609 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 6610 ! 6611 !-- Temporally comment out calculation of effective radiative temperature. 6612 !-- See below for more explanation. 6613 !-- (3) temperature 6614 !-- first we calculate an effective horizontal area to account for 6615 !-- the effect of vertical surfaces (which contributes to LW emission) 6616 !-- We simply use the ratio of the total LW to the incoming LW flux 6617 area_hor = pinlw / rad_lw_in_diff(nyn,nxl) 6618 t_rad_urb = ( ( pemitlw - pabslw + emissivity_urb * pinlw ) / & 6619 (emissivity_urb * sigma_sb * area_hor) )**0.25_wp 6633 6634 !- (2) emmsivity Eq. * in Krc et al. 2020 6635 !- emissivity_urb weighted average of surface and PC emissivity 6636 !- = absorbed LW in [surfaces + plant canopy] / pinlw 6637 IF ( pinsw /= 0.0_wp ) & 6638 emissivity_urb = (pabs_surf_lwdif + pabs_pc_lwdif) / pinlw 6639 6640 !- (3) temperature 6641 !- effective horizontal area to account for the effect of vertical 6642 !- surfaces, Eq. * in Krc et al. 2020 6643 area_norm = pinlw / rad_lw_in_diff(nyn,nxl) 6644 !- temperature, Eq. * in Krc et al. 2020 6645 t_rad_urb = ( (pemitlw - pabslw + emissivity_urb * pinlw) / & 6646 (emissivity_urb * sigma_sb * area_norm) )**0.25_wp 6620 6647 6621 6648 IF ( debug_output_timestep ) CALL debug_message( 'radiation_interaction', 'end' )
Note: See TracChangeset
for help on using the changeset viewer.
