Changeset 3378
- Timestamp:
- Oct 19, 2018 12:34:59 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
radiation_model_mod.f90 (modified) (17 diffs)
-
time_integration.f90 (modified) (2 diffs)
-
urban_surface_mod.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r3372 r3378 28 28 ! ----------------- 29 29 ! $Id$ 30 ! merge from radiation branch (r3362) into trunk 31 ! (moh.hefny): 32 ! - removed read/write_svf_on_init and read_dist_max_svf (not used anymore) 33 ! - bugfix nzut > nzpt in calculating maxboxes 34 ! 35 ! 3372 2018-10-18 14:03:19Z raasch 30 36 ! bugfix: kind type of 2nd argument of mpi_win_allocate changed, misplaced 31 37 ! __parallel directive … … 619 625 radiation_interactions = .FALSE., & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist) 620 626 surface_reflections = .TRUE., & !< flag to switch the calculation of radiation interaction between surfaces. 621 !< When it switched off, only the effect of buildings and trees shadow will627 !< When it switched off, only the effect of buildings and trees shadow 622 628 !< will be considered. However fewer SVFs are expected. 623 629 radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees … … 776 782 rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d) 777 783 rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 778 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m )779 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m )784 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m2) 785 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m2) 780 786 781 787 REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation … … 987 993 INTEGER(iwp), PARAMETER :: gasize = 100000 !< initial size of growing arrays 988 994 REAL(wp), PARAMETER :: grow_factor = 1.4_wp !< growth factor of growing arrays 989 REAL(wp) :: dist_max_svf = -9999.0 !< maximum distance to calculate the minimum svf to be considered. It is990 !< used to avoid very small SVFs resulting from too far surfaces with mutual visibility991 995 INTEGER(iwp) :: nsvfl !< number of svf for local processor 992 996 INTEGER(iwp) :: ncsfl !< no. of csf in local processor … … 1158 1162 skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,& 1159 1163 zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, & 1160 write_svf_on_init, read_svf_on_init, & 1161 nrefsteps, dist_max_svf, nsvfl, svf, & 1164 nrefsteps, nsvfl, svf, & 1162 1165 svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir, & 1163 1166 surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir, & … … 1167 1170 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 1168 1171 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 1169 nsurf_type, nzub, nzut, nz pt, nzu, pch, nsurf,&1172 nsurf_type, nzub, nzut, nzu, pch, nsurf, & 1170 1173 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 1171 1174 radiation_interactions, startwall, startland, endland, endwall, & … … 2910 2913 radiation_scheme, skip_time_do_radiation, & 2911 2914 sw_radiation, unscheduled_radiation_calls, & 2912 read_svf_on_init, write_svf_on_init, &2913 2915 max_raytracing_dist, min_irrf_value, & 2914 2916 nrefsteps, raytrace_mpi_rma, & 2915 dist_max_svf, &2916 2917 surface_reflections, svfnorm_report_thresh, & 2917 2918 radiation_interactions_on, & … … 2929 2930 max_raytracing_dist, min_irrf_value, & 2930 2931 nrefsteps, raytrace_mpi_rma, & 2931 dist_max_svf, &2932 2932 surface_reflections, svfnorm_report_thresh, & 2933 2933 radiation_interactions_on, & … … 4629 4629 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4630 4630 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 4631 ! REAL(wp) :: count_surfaces !< number of all surfaces in model domain4632 ! REAL(wp) :: count_surfaces_l !< number of all surfaces in sub-domain4633 ! REAL(wp) :: pt_surf_urb !< mean surface temperature of all surfaces in model domain, temporal work-around4634 ! REAL(wp) :: pt_surf_urb_l !< mean surface temperature of all surfaces in sub-domain, temporal work-around4635 4636 4631 REAL(wp), DIMENSION(0:nsurf_type) :: facearea 4637 4632 REAL(wp) :: pabsswl = 0.0_wp !< total absorbed SW radiation energy in local processor (W) … … 4686 4681 !-- precompute effective box depth with prototype Leaf Area Density 4687 4682 pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1 4688 CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift), &4683 CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift), & 4689 4684 60, prototype_lad, & 4690 4685 CSHIFT(ABS(sunorig), pc_box_dimshift), & 4691 4686 pc_box_area, pc_abs_frac) 4692 pc_box_area = pc_box_area * ABS(sunorig(pc_box_dimshift+1) / sunorig(1)) 4687 pc_box_area = pc_box_area * ABS(sunorig(pc_box_dimshift+1) & 4688 / sunorig(1)) 4693 4689 pc_abs_eff = LOG(1._wp - pc_abs_frac) / prototype_lad 4694 4690 ENDIF … … 4838 4834 j = surfl(iy, isurf) 4839 4835 i = surfl(ix, isurf) 4840 surfinswdir(isurf) = rad_sw_in_dir(j,i) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0) 4836 surfinswdir(isurf) = rad_sw_in_dir(j,i) * & 4837 costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0) 4841 4838 ENDDO 4842 4839 ! … … 5248 5245 t_rad_urb = ( (pemitlw - pabslw + emissivity_urb*pinlw) / & 5249 5246 (emissivity_urb*sigma_sb * area_hor) )**0.25_wp 5250 !5251 !-- It has been turned out that the effective radiative temperature is far5252 !-- too high during nighttime, resulting in unphysical radiative forcing5253 !-- with wrong signs. For the moment, as a work-around, compute the mean5254 !-- surface temperature from all surface elements, resulting in more5255 !-- physically meaningful radiative forcings.5256 ! pt_surf_urb_l = 0.0_wp5257 ! count_surfaces_l = 0.0_wp5258 ! DO m = 1, surf_lsm_h%ns5259 ! k = surf_lsm_h%k(m)5260 ! pt_surf_urb_l = pt_surf_urb_l + surf_lsm_h%pt_surface(m) &5261 ! * exner(k)5262 ! count_surfaces_l = count_surfaces_l + 1.0_wp5263 ! ENDDO5264 ! DO m = 1, surf_usm_h%ns5265 ! k = surf_usm_h%k(m)5266 ! pt_surf_urb_l = pt_surf_urb_l + surf_usm_h%pt_surface(m) &5267 ! * exner(k)5268 ! count_surfaces_l = count_surfaces_l + 1.0_wp5269 ! ENDDO5270 ! DO l = 0, 35271 ! DO m = 1, surf_lsm_v(l)%ns5272 ! k = surf_lsm_v(l)%k(m)5273 ! pt_surf_urb_l = pt_surf_urb_l + surf_lsm_v(l)%pt_surface(m) &5274 ! * exner(k)5275 ! count_surfaces_l = count_surfaces_l + 1.0_wp5276 ! ENDDO5277 ! DO m = 1, surf_usm_v(l)%ns5278 ! k = surf_usm_v(l)%k(m)5279 ! pt_surf_urb_l = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m) * exner(k)5280 ! count_surfaces_l = count_surfaces_l + 1.0_wp5281 ! ENDDO5282 ! ENDDO5283 !5284 ! pt_surf_urb = 0.0_wp5285 ! count_surfaces = 0.0_wp5286 !5287 ! #if defined( __parallel )5288 ! CALL MPI_ALLREDUCE( count_surfaces_l, count_surfaces, 1, MPI_REAL, &5289 ! MPI_SUM, comm2d, ierr)5290 ! CALL MPI_ALLREDUCE( pt_surf_urb_l, pt_surf_urb, 1, MPI_REAL, &5291 ! MPI_SUM, comm2d, ierr)5292 ! #else5293 ! count_surfaces_l = count_surfaces5294 ! pt_surf_urb_l = pt_surf_urb5295 ! #endif5296 !5297 ! t_rad_urb = pt_surf_urb / count_surfaces5298 5299 5247 5300 5248 CONTAINS … … 5469 5417 #endif 5470 5418 5471 5472 !INTEGER(iwp), DIMENSION(1:4,inorth_b:iwest_b) :: ijdb !< start and end of the local domain border coordinates (set in code)5473 !LOGICAL, DIMENSION(inorth_b:iwest_b) :: isborder !< is PE on the border of the domain in four corresponding directions5474 5475 5419 ! 5476 5420 !-- Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be … … 5566 5510 !-- check max_raytracing_dist relative to urban surface layer height 5567 5511 mrl = 2.0_wp * nzu * dz(1) 5512 !-- set max_raytracing_dist to double the urban surface layer height, if not set 5568 5513 IF ( max_raytracing_dist == -999.0_wp ) THEN 5569 5514 max_raytracing_dist = mrl 5515 ENDIF 5516 !-- check if max_raytracing_dist set too low (here we only warn the user. Other 5517 ! option is to correct the value again to double the urban surface layer height) 5518 IF ( max_raytracing_dist < mrl ) THEN 5519 WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist is set less than ', & 5520 'double the urban surface layer height, i.e. ', mrl 5521 CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0) 5570 5522 ENDIF 5571 5523 ! IF ( max_raytracing_dist <= mrl ) THEN … … 6303 6255 ENDDO 6304 6256 6305 6306 !--Advance itarget indices6257 ! 6258 !-- Advance itarget indices 6307 6259 itarg0 = itarg1 + 1 6308 6260 itarg1 = itarg1 + nzn … … 7467 7419 !-- 7468 7420 maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub, & 7469 nz pt - CEILING(origin(1)-.5_wp))) * nrays7421 nzut - CEILING(origin(1)-.5_wp))) * nrays 7470 7422 IF ( ncsfl + maxboxes > ncsfla ) THEN 7471 7423 !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) … … 8284 8236 !> It is called outside from usm_urban_surface_mod whenever the radiation fluxes 8285 8237 !> are needed. 8238 !> 8239 !> This routine is not used so far. However, it may serve as an interface for radiation 8240 !> fluxes of urban and land surfaces 8286 8241 !> 8287 8242 !> TODO: -
palm/trunk/SOURCE/time_integration.f90
r3347 r3378 25 25 ! ----------------- 26 26 ! $Id$ 27 ! merge from radiation branch (r3362) into trunk 28 ! (moh.hefny): 29 ! Bugfix in the if statement to call radiation_interaction 30 ! 31 ! 3347 2018-10-15 14:21:08Z suehring 27 32 ! - offline nesting separated from large-scale forcing module 28 33 ! - changes for synthetic turbulence generator … … 1149 1154 CALL cpu_log( log_point(50), 'radiation', 'stop' ) 1150 1155 1151 IF ( urban_surface .OR. land_surface .AND.&1156 IF ( ( urban_surface .OR. land_surface ) .AND. & 1152 1157 radiation_interactions ) THEN 1153 1158 CALL cpu_log( log_point(75), 'radiation_interaction', 'start' ) -
palm/trunk/SOURCE/urban_surface_mod.f90
r3371 r3378 28 28 ! ----------------- 29 29 ! $Id$ 30 ! merge from radiation branch (r3362) into trunk 31 ! (moh.hefny): 32 ! - check the requested output variables if they are correct 33 ! - added unscheduled_radiation_calls switch to control force_radiation_call 34 ! - minor formate changes 35 ! 36 ! 3371 2018-10-18 13:40:12Z knoop 30 37 ! Set flag indicating that albedo at urban surfaces is already initialized 31 38 ! … … 399 406 idcsf, ndcsf, kdcsf, pct, & 400 407 startland, endland, startwall, endwall, skyvf, skyvft, nzub, & 401 nzut, n zpt, npcbl, pcbl408 nzut, npcbl, pcbl, unscheduled_radiation_calls 402 409 403 410 USE statistics, & … … 2399 2406 IMPLICIT NONE 2400 2407 2401 CHARACTER (len=*),INTENT(IN) :: variable !: 2402 CHARACTER (len=*),INTENT(OUT) :: unit !: 2403 2404 CHARACTER (len=varnamelength) :: var 2405 2408 CHARACTER(LEN=*),INTENT(IN) :: variable !: 2409 CHARACTER(LEN=*),INTENT(OUT) :: unit !: 2410 2411 INTEGER(iwp) :: i,j !< index 2412 CHARACTER(LEN=varnamelength) :: var 2413 INTEGER(iwp), PARAMETER :: nl1 = 32 !< number of directional usm variables 2414 CHARACTER(LEN=varnamelength), DIMENSION(nl1) :: varlist1 = & !< list of directional usm variables 2415 (/'usm_rad_net','usm_rad_insw','usm_rad_inlw','usm_rad_inswdir', & 2416 'usm_rad_inswdif','usm_rad_inswref','usm_rad_inlwdif','usm_wshf', & 2417 'usm_rad_inlwref','usm_rad_outsw','usm_rad_outlw','usm_rad_hf', & 2418 'usm_rad_ressw','usm_rad_reslw','usm_wghf','usm_wghf_window', & 2419 'usm_wghf_green','usm_iwghf','usm_iwghf_window','usm_surfz', & 2420 'usm_surfwintrans','usm_surfcat','usm_surfalb','usm_surfemis', & 2421 'usm_t_surf','usm_t_surf_window','usm_t_surf_green','usm_t_green', & 2422 'usm_t_surf_10cm','usm_t_wall','usm_t_window','usm_t_green'/) 2423 2424 INTEGER(iwp), PARAMETER :: nl2 = 9 !< number of other variables 2425 CHARACTER(LEN=varnamelength), DIMENSION(nl2) :: varlist2 = & !< list of other usm variables 2426 (/'usm_skyvf','usm_skyvft','usm_svf','usm_dif','usm_rad_pc_inlw', & 2427 'usm_rad_pc_insw','usm_rad_pc_inswdir','usm_rad_pc_inswdif', & 2428 'usm_rad_pc_inswref'/) 2429 2430 INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions 2431 CHARACTER(LEN=6), DIMENSION(nd), PARAMETER :: dirname = & !< direction names 2432 (/'_roof ','_south','_north','_west ','_east '/) 2433 LOGICAL :: lfound !< flag if the variable is found 2434 2435 2436 lfound = .FALSE. 2406 2437 var = TRIM(variable) 2438 2439 !-- check if variable exists 2440 ! directional variables 2441 DO i = 1, nl1 2442 DO j = 1, nd 2443 IF ( TRIM(var) == TRIM(varlist1(i))//TRIM(dirname(j)) ) THEN 2444 lfound = .TRUE. 2445 EXIT 2446 ENDIF 2447 IF ( lfound ) EXIT 2448 ENDDO 2449 ENDDO 2450 IF ( lfound ) GOTO 10 2451 ! other variables 2452 DO i = 1, nl2 2453 IF ( TRIM(var) == TRIM(varlist2(i)) ) THEN 2454 lfound = .TRUE. 2455 EXIT 2456 ENDIF 2457 ENDDO 2458 IF ( .NOT. lfound ) THEN 2459 unit = 'illegal' 2460 RETURN 2461 ENDIF 2462 10 CONTINUE 2463 2464 !-- check if variable exists 2407 2465 IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. & 2408 2466 var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. & … … 2593 2651 2594 2652 CASE ( 'usm_surfz' ) 2595 !-- array of lw radiation falling to local surface after i-th reflection2653 !-- array of surface height (z) 2596 2654 IF ( idsint == iup_u ) THEN 2597 2655 DO m = 1, surf_usm_h%ns … … 5332 5390 ENDIF 5333 5391 ENDDO 5334 !!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!5335 ! t_window_v_p(l)%t = t_wall_v_p(l)%t5336 ! surf_usm_v(l)%tt_window_m = surf_usm_v(l)%tt_wall_m5337 ! t_green_v_p(l)%t = t_wall_v_p(l)%t5338 ! surf_usm_v(l)%tt_green_m = surf_usm_v(l)%tt_wall_m5339 !!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!5340 5392 ENDDO 5341 5393 … … 6914 6966 INTEGER(iwp) :: wtc 6915 6967 REAL(wp), DIMENSION(n_surface_params) :: wtp 6916 6917 6968 LOGICAL :: ascii_file = .FALSE. 6918 6919 6969 INTEGER(iwp), DIMENSION(0:17, nysg:nyng, nxlg:nxrg) :: usm_par 6920 6970 REAL(wp), DIMENSION(1:14, nysg:nyng, nxlg:nxrg) :: usm_val … … 7742 7792 !-- in case of fast changes in the skin temperature, it is required to 7743 7793 !-- update the radiative fluxes in order to keep the solution stable 7744 IF ( ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) .OR.&7794 IF ( ( ( ABS( t_surf_h_p(m) - t_surf_h(m) ) > 1.0_wp ) .OR. & 7745 7795 ( ABS( t_surf_green_h_p(m) - t_surf_green_h(m) ) > 1.0_wp ) .OR. & 7746 ( ABS( t_surf_window_h_p(m) - t_surf_window_h(m) ) > 1.0_wp ) ) THEN 7796 ( ABS( t_surf_window_h_p(m) - t_surf_window_h(m) ) > 1.0_wp ) ) & 7797 .AND. unscheduled_radiation_calls ) THEN 7747 7798 force_radiation_call_l = .TRUE. 7748 7799 ENDIF … … 7897 7948 7898 7949 !-- add RK3 term 7899 t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) * &7950 t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) * & 7900 7951 surf_usm_v(l)%tt_surface_m(m) 7901 7952 t_surf_window_v_p(l)%t(m) = t_surf_window_v_p(l)%t(m) + dt_3d * tsc(3) * & … … 7907 7958 !-- store also vpt_surface, which is, due to the lack of moisture on roofs simply 7908 7959 !-- assumed to be the surface temperature. 7909 surf_usm_v(l)%pt_surface(m) = ( surf_usm_v(l)%frac(ind_veg_wall,m) * t_surf_v_p(l)%t(m) &7910 + surf_usm_v(l)%frac(ind_wat_win,m) * t_surf_window_v_p(l)%t(m) &7960 surf_usm_v(l)%pt_surface(m) = ( surf_usm_v(l)%frac(ind_veg_wall,m) * t_surf_v_p(l)%t(m) & 7961 + surf_usm_v(l)%frac(ind_wat_win,m) * t_surf_window_v_p(l)%t(m) & 7911 7962 + surf_usm_v(l)%frac(ind_pav_green,m) * t_surf_green_v_p(l)%t(m) ) & 7912 7963 / exner(k) … … 7920 7971 stend_window = ( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) - dt_3d * tsc(3) *& 7921 7972 surf_usm_v(l)%tt_surface_window_m(m) ) / ( dt_3d * tsc(2) ) 7922 stend_green = ( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) - dt_3d * tsc(3) * &7973 stend_green = ( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) - dt_3d * tsc(3) * & 7923 7974 surf_usm_v(l)%tt_surface_green_m(m) ) / ( dt_3d * tsc(2) ) 7924 7975 … … 7929 7980 surf_usm_v(l)%tt_surface_window_m(m) = stend_window 7930 7981 surf_usm_v(l)%tt_surface_green_m(m) = stend_green 7931 ELSEIF ( intermediate_timestep_count < &7982 ELSEIF ( intermediate_timestep_count < & 7932 7983 intermediate_timestep_count_max ) THEN 7933 surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend + &7984 surf_usm_v(l)%tt_surface_m(m) = -9.5625_wp * stend + & 7934 7985 5.3125_wp * surf_usm_v(l)%tt_surface_m(m) 7935 surf_usm_v(l)%tt_surface_green_m(m) = -9.5625_wp * stend_green + 7986 surf_usm_v(l)%tt_surface_green_m(m) = -9.5625_wp * stend_green + & 7936 7987 5.3125_wp * surf_usm_v(l)%tt_surface_green_m(m) 7937 surf_usm_v(l)%tt_surface_window_m(m) = -9.5625_wp * stend_window + 7988 surf_usm_v(l)%tt_surface_window_m(m) = -9.5625_wp * stend_window + & 7938 7989 5.3125_wp * surf_usm_v(l)%tt_surface_window_m(m) 7939 7990 ENDIF … … 7943 7994 !-- update the radiative fluxes in order to keep the solution stable 7944 7995 7945 IF ( ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) .OR.&7996 IF ( ( ( ABS( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) ) > 1.0_wp ) .OR. & 7946 7997 ( ABS( t_surf_green_v_p(l)%t(m) - t_surf_green_v(l)%t(m) ) > 1.0_wp ) .OR. & 7947 ( ABS( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) ) > 1.0_wp ) ) THEN 7998 ( ABS( t_surf_window_v_p(l)%t(m) - t_surf_window_v(l)%t(m) ) > 1.0_wp ) ) & 7999 .AND. unscheduled_radiation_calls ) THEN 7948 8000 force_radiation_call_l = .TRUE. 7949 8001 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
