Changeset 2977 for palm/trunk
- Timestamp:
- Apr 17, 2018 10:27:57 AM (7 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 7 edited
-
init_3d_model.f90 (modified) (5 diffs)
-
palm.f90 (modified) (3 diffs)
-
plant_canopy_model_mod.f90 (modified) (9 diffs)
-
radiation_model_mod.f90 (modified) (23 diffs)
-
surface_mod.f90 (modified) (6 diffs)
-
time_integration.f90 (modified) (2 diffs)
-
urban_surface_mod.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_3d_model.f90
r2938 r2977 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implement changes from branch radiation (r2948-2971) with minor modifications 28 ! (moh.hefny): 29 ! - set radiation_interactions according to the existence of urban/land vertical 30 ! surfaces and trees to activiate RTM 31 ! - set average_radiation to TRUE if RTM is activiated 32 ! 33 ! 2938 2018-03-27 15:52:42Z suehring 27 34 ! - Revise Inifor initialization for geostrophic wind components 28 35 ! - Initialize synthetic turbulence generator in case of Inifor initialization … … 506 513 507 514 USE radiation_model_mod, & 508 ONLY: radiation_init, radiation, radiation_control, radiation_scheme, & 515 ONLY: average_radiation, & 516 radiation_init, radiation, radiation_control, radiation_scheme, & 509 517 radiation_calc_svf, radiation_write_svf, & 510 518 radiation_interaction, radiation_interactions, & 511 519 radiation_interaction_init, radiation_read_svf, & 512 radiation_presimulate_solar_pos 520 radiation_presimulate_solar_pos, radiation_interactions_on 513 521 514 522 USE random_function_mod … … 534 542 USE surface_mod, & 535 543 ONLY : init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h, & 536 surf_usm_h, get_topography_top_index_ji 544 surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist 537 545 538 546 USE transpose_indices … … 2357 2365 !-- If required, set chemical emissions 2358 2366 !-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model) 2359 IF ( air_chemistry ) THEN2367 IF ( air_chemistry ) THEN 2360 2368 CALL chem_emissions 2361 2369 ENDIF … … 2364 2372 !-- Initialize radiation processes 2365 2373 IF ( radiation ) THEN 2374 ! 2375 !-- Activate radiation_interactions according to the existence of vertical surfaces and/or trees. 2376 !-- The namelist parameter radiation_interactions_on can override this behavior. 2377 !-- (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 2378 !-- init_surface_arrays.) 2379 IF ( radiation_interactions_on ) THEN 2380 IF ( vertical_surfaces_exist .OR. plant_canopy ) THEN 2381 radiation_interactions = .TRUE. 2382 average_radiation = .TRUE. 2383 ELSE 2384 radiation_interactions_on = .FALSE. !< reset namelist parameter: no interactions 2385 !< calculations necessary in case of flat surface 2386 ENDIF 2387 ELSEIF ( vertical_surfaces_exist .OR. plant_canopy ) THEN 2388 message_string = 'radiation_interactions_on is set to .FALSE. although ' // & 2389 'vertical surfaces and/or trees exist. The model will run ' // & 2390 'without RTM (no shadows, no radiation reflections)' 2391 CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 ) 2392 ENDIF 2366 2393 ! 2367 2394 !-- If required, initialize radiation interactions between surfaces -
palm/trunk/SOURCE/palm.f90
r2951 r2977 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Deduct spinup_time from RUN_CONTROL output of main 3d run 28 ! (use time_since_reference_point instead of simulated_time) 29 ! 30 ! 2951 2018-04-06 09:05:08Z kanani 27 31 ! Add log_point_s for pmci_init 28 32 ! … … 240 244 nest_domain, neutral, nudging, passive_scalar, runnr, & 241 245 simulated_time, simulated_time_chr, spinup, & 246 time_since_reference_point, & 242 247 user_interface_current_revision, & 243 248 user_interface_required_revision, version, wall_heatflux, & … … 464 469 ! 465 470 !-- Set start time in format hh:mm:ss 466 simulated_time_chr = time_to_string( simulated_time)471 simulated_time_chr = time_to_string( time_since_reference_point ) 467 472 468 473 ! -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r2932 r2977 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- 21 ! ------------------ 22 22 ! 23 23 ! … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implement changes from branch radiation (r2948-2971) with minor modifications, 28 ! plus some formatting. 29 ! (moh.hefny): 30 ! Add plant canopy type to account for changes in LAD (based on the changes 31 ! done by Resler & Pavel) and correct the error message to PALM Standard. 32 ! 33 ! 2932 2018-03-26 09:39:22Z maronga 27 34 ! renamed canopy_par to plant_canopy_parameters 28 35 ! … … 213 220 REAL(wp) :: lad_vertical_gradient(10) = 0.0_wp !< lad gradient 214 221 REAL(wp) :: lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m) 222 223 REAL(wp) :: lad_type_coef(0:10) = 1.0_wp !< multiplicative coeficients for particular types 224 !< of plant canopy (e.g. deciduous tree during winter) 215 225 216 226 REAL(wp), DIMENSION(:), ALLOCATABLE :: lad !< leaf area density … … 962 972 alpha_lad, beta_lad, canopy_drag_coeff, & 963 973 canopy_mode, cthf, & 964 lad_surface, 974 lad_surface, lad_type_coef, & 965 975 lad_vertical_gradient, & 966 976 lad_vertical_gradient_level, & … … 971 981 NAMELIST /canopy_par/ alpha_lad, beta_lad, canopy_drag_coeff, & 972 982 canopy_mode, cthf, & 973 lad_surface, 983 lad_surface, lad_type_coef, & 974 984 lad_vertical_gradient, & 975 985 lad_vertical_gradient_level, & … … 1034 1044 !> 1035 1045 !> num_levels 1036 !> dtype,x,y, value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)1037 !> dtype,x,y, value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)1038 !> dtype,x,y, value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)1046 !> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1) 1047 !> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1) 1048 !> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1) 1039 1049 !> ... 1040 1050 !> … … 1060 1070 1061 1071 INTEGER(iwp) :: dtype !< type of input data (1=lad) 1072 INTEGER(iwp) :: pctype !< type of plant canopy (deciduous,non-deciduous,...) 1062 1073 INTEGER(iwp) :: i, j !< running index 1063 1074 INTEGER(iwp) :: nzp !< number of vertical layers of plant canopy … … 1073 1084 ! 1074 1085 !-- Open and read plant canopy input data 1075 OPEN(152, file='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ), &1076 access='SEQUENTIAL', action='READ', status='OLD', &1077 form='FORMATTED', err=515)1078 READ(152, *, err=516, end=517)nzp !< read first line = number of vertical layers1086 OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ), & 1087 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD', & 1088 FORM='FORMATTED', ERR=515) 1089 READ(152, *, ERR=516, END=517) nzp !< read first line = number of vertical layers 1079 1090 1080 ALLOCATE( col(0:nzp-1))1091 ALLOCATE( col(0:nzp-1) ) 1081 1092 1082 1093 DO 1083 READ(152, *, err=516, end=517) dtype, i, j, col(:) 1084 IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE 1085 1086 SELECT CASE (dtype) 1087 CASE( 1 ) !< leaf area density 1088 ! 1089 !-- This is just the pure canopy layer assumed to be grounded to 1090 !-- a flat domain surface. At locations where plant canopy sits 1091 !-- on top of any kind of topography, the vertical plant column 1092 !-- must be "lifted", which is done in SUBROUTINE pcm_tendency. 1093 lad_s(0:nzp-1, j, i) = col(0:nzp-1) 1094 1095 CASE DEFAULT 1096 WRITE(message_string, '(a,i2,a)') & 1097 'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"' 1098 CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 ) 1099 END SELECT 1094 READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:) 1095 IF ( i < nxlg .OR. i > nxrg .OR. j < nysg .OR. j > nyng ) CYCLE 1096 1097 SELECT CASE (dtype) 1098 CASE( 1 ) !< leaf area density 1099 ! 1100 !-- This is just the pure canopy layer assumed to be grounded to 1101 !-- a flat domain surface. At locations where plant canopy sits 1102 !-- on top of any kind of topography, the vertical plant column 1103 !-- must be "lifted", which is done in SUBROUTINE pcm_tendency. 1104 IF ( pctype < 0 .OR. pctype > 10 ) THEN !< incorrect plant canopy type 1105 WRITE( message_string, * ) 'Incorrect type of plant canopy. ' // & 1106 'Allowed values 0 <= pctype <= 10, ' // & 1107 'but pctype is ', pctype 1108 CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 ) 1109 ENDIF 1110 lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype) 1111 1112 CASE DEFAULT 1113 WRITE(message_string, '(a,i2,a)') & 1114 'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"' 1115 CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 ) 1116 END SELECT 1100 1117 ENDDO 1101 1118 … … 1107 1124 1108 1125 517 CLOSE(152) 1109 DEALLOCATE( col)1126 DEALLOCATE( col ) 1110 1127 1111 1128 CALL exchange_horiz( lad_s, nbgp ) -
palm/trunk/SOURCE/radiation_model_mod.f90
r2967 r2977 22 22 ! 23 23 ! Current revisions: 24 ! ----------------- 24 ! ------------------ 25 25 ! 26 26 ! … … 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Implement changes from branch radiation (r2948-2971) with minor modifications, 31 ! plus some formatting. 32 ! (moh.hefny): 33 ! - replaced plant_canopy by npcbl to check tree existence to avoid weird 34 ! allocation of related arrays (after domain decomposition some domains 35 ! contains no trees although plant_canopy (global parameter) is still TRUE). 36 ! - added a namelist parameter to force RTM settings 37 ! - enabled the option to switch radiation reflections off 38 ! - renamed surf_reflections to surface_reflections 39 ! - removed average_radiation flag from the namelist (now it is implicitly set 40 ! in init_3d_model according to RTM) 41 ! - edited read and write sky view factors and CSF routines to account for 42 ! the sub-domains which may not contain any of them 43 ! 44 ! 2967 2018-04-13 11:22:08Z raasch 30 45 ! bugfix: missing parallel cpp-directives added 31 46 ! … … 432 447 sun_direction = .FALSE., & !< flag parameter indicating whether solar direction shall be calculated 433 448 average_radiation = .FALSE., & !< flag to set the calculation of radiation averaging for the domain 434 radiation_interactions = .FALSE., & !< flag to control if radiation interactions via sky-view factors shall be considered435 surf _reflections = .TRUE.!< flag to switch the calculation of radiation interaction between surfaces.449 radiation_interactions = .FALSE., & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist) 450 surface_reflections = .TRUE., & !< flag to switch the calculation of radiation interaction between surfaces. 436 451 !< When it switched off, only the effect of buildings and trees shadow will 437 452 !< will be considered. However fewer SVFs are expected. 453 radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees 438 454 439 455 … … 680 696 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy 681 697 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy 682 INTEGER(iwp) :: npcbl 698 INTEGER(iwp) :: npcbl = 0 !< number of the plant canopy gridboxes in local processor 683 699 INTEGER(wp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i] 684 700 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box … … 955 971 surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart, & 956 972 surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw, & 957 pcbl, npcbl,iup_u, inorth_u, isouth_u, ieast_u, iwest_u, &973 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 958 974 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 959 975 nsurf_type, nzub, nzut, nzu, pch, nsurf, & … … 961 977 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 962 978 radiation_interactions, startwall, startland, endland, endwall, & 963 skyvf, skyvft 979 skyvf, skyvft, radiation_interactions_on, average_radiation 964 980 965 981 #if defined ( __rrtmg ) … … 1309 1325 CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 ) 1310 1326 ENDIF 1311 ENDIF1312 1313 !1314 !-- Radiation interactions1315 IF ( nrefsteps < 1 .AND. radiation_interactions ) THEN1316 message_string = 'nrefsteps must be > 0 when using LSM/USM to' // &1317 'account for surface outgoing SW flux.' // &1318 'You may set surf_reflections = .FALSE. to ' // &1319 'diable surface reflections instead.'1320 CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )1321 1327 ENDIF 1322 1328 … … 2634 2640 nrefsteps, mrt_factors, rma_lad_raytrace, & 2635 2641 dist_max_svf, & 2636 average_radiation,&2637 surf_reflections, svfnorm_report_thresh2642 surface_reflections, svfnorm_report_thresh, & 2643 radiation_interactions_on 2638 2644 2639 2645 NAMELIST /radiation_parameters/ albedo, albedo_type, albedo_lw_dir, & … … 2647 2653 nrefsteps, mrt_factors, rma_lad_raytrace, & 2648 2654 dist_max_svf, & 2649 average_radiation,&2650 surf_reflections, svfnorm_report_thresh2655 surface_reflections, svfnorm_report_thresh, & 2656 radiation_interactions_on 2651 2657 2652 2658 line = ' ' … … 2693 2699 radiation = .TRUE. 2694 2700 2701 IF ( .NOT. radiation_interactions_on .AND. surface_reflections ) THEN 2702 message_string = 'surface_reflections is allowed only when ' // & 2703 'radiation_interactions_on is set to TRUE' 2704 CALL message( 'radiation_parin', 'PA0487',1, 2, 0, 6, 0 ) 2705 ENDIF 2706 2695 2707 12 CONTINUE 2696 2697 2698 2699 !-- Set radiation_interactions flag according to urban_ and land_surface flag2700 IF ( urban_surface .OR. land_surface ) radiation_interactions = .TRUE.2701 2708 2702 2709 END SUBROUTINE radiation_parin … … 4354 4361 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) 4355 4362 4356 IF ( plant_canopy) THEN4363 IF ( npcbl > 0 ) THEN 4357 4364 !-- precompute effective box depth with prototype Leaf Area Density 4358 4365 pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1 … … 4465 4472 #endif 4466 4473 4467 DO isvf = 1, nsvfl 4468 isurf = svfsurf(1, isvf) 4469 k = surfl(iz, isurf) 4470 j = surfl(iy, isurf) 4471 i = surfl(ix, isurf) 4472 isurfsrc = svfsurf(2, isvf) 4473 4474 !-- for surface-to-surface factors we calculate thermal radiation in 1st pass 4475 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4476 ENDDO 4474 IF ( surface_reflections) THEN 4475 DO isvf = 1, nsvfl 4476 isurf = svfsurf(1, isvf) 4477 k = surfl(iz, isurf) 4478 j = surfl(iy, isurf) 4479 i = surfl(ix, isurf) 4480 isurfsrc = svfsurf(2, isvf) 4481 ! 4482 !-- For surface-to-surface factors we calculate thermal radiation in 1st pass 4483 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4484 ENDDO 4485 ENDIF 4477 4486 4478 4487 !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors … … 4494 4503 ENDIF 4495 4504 4496 IF ( plant_canopy) THEN4505 IF ( npcbl > 0 ) THEN 4497 4506 4498 4507 pcbinswdir(:) = 0._wp … … 4536 4545 ! surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw 4537 4546 4547 IF ( .NOT. surface_reflections ) THEN 4548 ! 4549 !-- Set nrefsteps to 0 to disable reflections 4550 nrefsteps = 0 4551 surfoutsl = albedo_surf * surfins 4552 surfoutll = (1._wp - emiss_surf) * surfinl 4553 surfoutsw = surfoutsw + surfoutsl 4554 surfoutlw = surfoutlw + surfoutll 4555 ENDIF 4556 4538 4557 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4539 4558 !-- Next passes - reflections … … 4585 4604 4586 4605 !-- push heat flux absorbed by plant canopy to respective 3D arrays 4587 IF ( plant_canopy) THEN4606 IF ( npcbl > 0 ) THEN 4588 4607 pc_heating_rate(:,:,:) = 0._wp 4589 4608 DO ipcgb = 1, npcbl … … 4733 4752 !-- domain when using average_radiation in the respective radiation model 4734 4753 4735 IF ( average_radiation ) THEN 4736 4737 !-- precalculate face areas for all face directions using normal vector 4738 DO d = 0, nsurf_type 4739 facearea(d) = 1._wp 4740 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 4741 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 4742 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 4743 ENDDO 4744 ! 4745 !-- absorbed/received SW & LW and emitted LW energy of all physical 4746 !-- surfaces (land and urban) in local processor 4747 pinswl = 0._wp 4748 pinlwl = 0._wp 4749 pabsswl = 0._wp 4750 pabslwl = 0._wp 4751 pemitlwl = 0._wp 4752 emiss_sum_surfl = 0._wp 4753 area_surfl = 0._wp 4754 DO i = 1, nsurfl 4755 d = surfl(id, i) 4756 !-- received SW & LW 4757 pinswl = pinswl + surfinsw(i) * facearea(d) 4758 pinlwl = pinlwl + surfinlw(i) * facearea(d) 4759 !-- absorbed SW & LW 4760 pabsswl = pabsswl + (1._wp - albedo_surf(i)) * & 4761 surfinsw(i) * facearea(d) 4762 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d) 4763 !-- emitted LW 4764 pemitlwl = pemitlwl + surfoutlw(i) * facearea(d) 4765 !-- emissivity and area sum 4766 emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d) 4767 area_surfl = area_surfl + facearea(d) 4768 END DO 4769 ! 4770 !-- add the absorbed SW energy by plant canopy 4771 IF ( plant_canopy ) THEN 4772 pabsswl = pabsswl + SUM(pcbinsw) 4773 pabslwl = pabslwl + SUM(pcbinlw) 4774 ENDIF 4775 ! 4776 !-- gather all rad flux energy in all processors 4754 !-- Precalculate face areas for all face directions using normal vector 4755 DO d = 0, nsurf_type 4756 facearea(d) = 1._wp 4757 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 4758 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 4759 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 4760 ENDDO 4761 ! 4762 !-- absorbed/received SW & LW and emitted LW energy of all physical 4763 !-- surfaces (land and urban) in local processor 4764 pinswl = 0._wp 4765 pinlwl = 0._wp 4766 pabsswl = 0._wp 4767 pabslwl = 0._wp 4768 pemitlwl = 0._wp 4769 emiss_sum_surfl = 0._wp 4770 area_surfl = 0._wp 4771 DO i = 1, nsurfl 4772 d = surfl(id, i) 4773 !-- received SW & LW 4774 pinswl = pinswl + surfinsw(i) * facearea(d) 4775 pinlwl = pinlwl + surfinlw(i) * facearea(d) 4776 !-- absorbed SW & LW 4777 pabsswl = pabsswl + (1._wp - albedo_surf(i)) * & 4778 surfinsw(i) * facearea(d) 4779 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d) 4780 !-- emitted LW 4781 pemitlwl = pemitlwl + surfoutlw(i) * facearea(d) 4782 !-- emissivity and area sum 4783 emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d) 4784 area_surfl = area_surfl + facearea(d) 4785 END DO 4786 ! 4787 !-- add the absorbed SW energy by plant canopy 4788 IF ( npcbl > 0 ) THEN 4789 pabsswl = pabsswl + SUM(pcbinsw) 4790 pabslwl = pabslwl + SUM(pcbinlw) 4791 ENDIF 4792 ! 4793 !-- gather all rad flux energy in all processors 4777 4794 #if defined( __parallel ) 4778 4779 4780 4781 4782 4783 4784 4795 CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 4796 CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 4797 CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4798 CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4799 CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4800 CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4801 CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4785 4802 #else 4786 4787 4788 4789 4790 4791 4792 4803 pinsw = pinswl 4804 pinlw = pinlwl 4805 pabssw = pabsswl 4806 pabslwl = pabslw 4807 pemitlwl = pemitlw 4808 emiss_sum_surf = emiss_sum_surfl 4809 area_surf = area_surfl 4793 4810 #endif 4794 4811 4795 !-- (1) albedo 4796 IF ( pinsw /= 0.0_wp ) albedo_urb = 1._wp - pabssw / pinsw 4797 4798 !-- (2) average emmsivity 4799 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 4800 4801 !-- (3) temperature 4802 t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp 4803 4804 ENDIF 4812 !-- (1) albedo 4813 IF ( pinsw /= 0.0_wp ) albedo_urb = 1._wp - pabssw / pinsw 4814 4815 !-- (2) average emmsivity 4816 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 4817 4818 !-- (3) temperature 4819 t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp 4820 4805 4821 4806 4822 CONTAINS … … 5071 5087 5072 5088 !-- fill gridpcbl and pcbl 5073 IF ( plant_canopy) THEN5089 IF ( npcbl > 0 ) THEN 5074 5090 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 5075 5091 ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) ) … … 5593 5609 5594 5610 DEALLOCATE ( zdirs, zbdry, vffrac, ztransp ) 5595 5596 DO isurfs = 1, nsurf 5597 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 5598 surfl(iz, isurflt), surfl(id, isurflt), & 5599 surf(ix, isurfs), surf(iy, isurfs), & 5600 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5601 CYCLE 5602 ENDIF 5611 ! 5612 !-- Following calculations only required for surface_reflections 5613 IF ( surface_reflections ) THEN 5614 5615 DO isurfs = 1, nsurf 5616 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 5617 surfl(iz, isurflt), surfl(id, isurflt), & 5618 surf(ix, isurfs), surf(iy, isurfs), & 5619 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5620 CYCLE 5621 ENDIF 5603 5622 5604 sd = surf(id, isurfs)5605 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &5606 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &5607 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)5608 5609 !-- unit vector source -> target5610 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)5611 sqdist = SUM(uv(:)**2)5612 uv = uv / SQRT(sqdist)5613 5614 !-- reject raytracing above max distance5615 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN5616 ray_skip_maxdist = ray_skip_maxdist + 15617 CYCLE5618 ENDIF5619 5620 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area5621 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction5622 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction5623 / (pi * sqdist) & ! square of distance between centers5624 * facearea(sd)5625 5626 !-- reject raytracing for potentially too small view factor values5627 IF ( rirrf < min_irrf_value ) THEN5628 ray_skip_minval = ray_skip_minval + 15629 CYCLE5630 ENDIF5631 5632 !-- raytrace + process plant canopy sinks within5633 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &5634 visible, transparency, win_lad)5635 5636 IF ( .NOT. visible ) CYCLE5637 ! rsvf = rirrf * transparency5638 5639 !-- write to the svf array5640 nsvfl = nsvfl + 15641 !-- check dimmension of asvf array and enlarge it if needed5642 IF ( nsvfla < nsvfl ) THEN5643 k = nsvfla * 25644 IF ( msvf == 0 ) THEN5623 sd = surf(id, isurfs) 5624 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 5625 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 5626 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 5627 5628 !-- unit vector source -> target 5629 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5630 sqdist = SUM(uv(:)**2) 5631 uv = uv / SQRT(sqdist) 5632 5633 !-- reject raytracing above max distance 5634 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN 5635 ray_skip_maxdist = ray_skip_maxdist + 1 5636 CYCLE 5637 ENDIF 5638 5639 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 5640 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 5641 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 5642 / (pi * sqdist) & ! square of distance between centers 5643 * facearea(sd) 5644 5645 !-- reject raytracing for potentially too small view factor values 5646 IF ( rirrf < min_irrf_value ) THEN 5647 ray_skip_minval = ray_skip_minval + 1 5648 CYCLE 5649 ENDIF 5650 5651 !-- raytrace + process plant canopy sinks within 5652 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., & 5653 visible, transparency, win_lad) 5654 5655 IF ( .NOT. visible ) CYCLE 5656 ! rsvf = rirrf * transparency 5657 5658 !-- write to the svf array 5659 nsvfl = nsvfl + 1 5660 !-- check dimmension of asvf array and enlarge it if needed 5661 IF ( nsvfla < nsvfl ) THEN 5662 k = nsvfla * 2 5663 IF ( msvf == 0 ) THEN 5645 5664 msvf = 1 5646 5665 ALLOCATE( asvf1(k) ) … … 5648 5667 asvf1(1:nsvfla) = asvf2 5649 5668 DEALLOCATE( asvf2 ) 5650 ELSE5669 ELSE 5651 5670 msvf = 0 5652 5671 ALLOCATE( asvf2(k) ) … … 5654 5673 asvf2(1:nsvfla) = asvf1 5655 5674 DEALLOCATE( asvf1 ) 5656 ENDIF 5657 5658 ! WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 5659 ! CALL radiation_write_debug_log( msg ) 5660 5661 nsvfla = k 5662 ENDIF 5663 !-- write svf values into the array 5664 asvf(nsvfl)%isurflt = isurflt 5665 asvf(nsvfl)%isurfs = isurfs 5666 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency 5667 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor 5668 ENDDO 5675 ENDIF 5676 5677 ! WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 5678 ! CALL radiation_write_debug_log( msg ) 5679 5680 nsvfla = k 5681 ENDIF 5682 !-- write svf values into the array 5683 asvf(nsvfl)%isurflt = isurflt 5684 asvf(nsvfl)%isurfs = isurfs 5685 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency 5686 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor 5687 ENDDO 5688 ENDIF 5669 5689 ENDDO 5670 5690 … … 6733 6753 !-- read nsvfl, ncsfl 6734 6754 READ ( fsvf ) nsvfl, ncsfl, nsurfl_from_file 6735 IF ( nsvfl < =0 .OR. ncsfl < 0 ) THEN6755 IF ( nsvfl < 0 .OR. ncsfl < 0 ) THEN 6736 6756 WRITE( message_string, * ) 'Wrong number of SVF or CSF' 6737 6757 CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 ) … … 6747 6767 CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 ) 6748 6768 ENDIF 6749 6750 IF ( .NOT. ALLOCATED( skyvf ) ) ALLOCATE( skyvf(nsurfl) ) 6751 IF ( .NOT. ALLOCATED( skyvft ) ) ALLOCATE( skyvft(nsurfl) ) 6752 IF ( .NOT. ALLOCATED( svf ) ) ALLOCATE( svf(ndsvf,nsvfl) ) 6753 IF ( .NOT. ALLOCATED( svfsurf ) ) ALLOCATE( svfsurf(idsvf,nsvfl) ) 6754 READ(fsvf) skyvf 6755 READ(fsvf) skyvft 6756 READ(fsvf) svf 6757 READ(fsvf) svfsurf 6758 IF ( plant_canopy ) THEN 6759 IF ( .NOT. ALLOCATED( csf ) ) ALLOCATE( csf(ndcsf,ncsfl) ) 6760 IF ( .NOT. ALLOCATED( csfsurf ) ) ALLOCATE( csfsurf(idcsf,ncsfl) ) 6761 READ(fsvf) csf 6762 READ(fsvf) csfsurf 6769 IF ( nsurfl > 0 ) THEN 6770 IF ( .NOT. ALLOCATED( skyvf ) ) ALLOCATE( skyvf(nsurfl) ) 6771 IF ( .NOT. ALLOCATED( skyvft ) ) ALLOCATE( skyvft(nsurfl) ) 6772 READ(fsvf) skyvf 6773 READ(fsvf) skyvft 6774 ENDIF 6775 6776 IF ( nsvfl > 0 ) THEN 6777 IF ( .NOT. ALLOCATED( svf ) ) ALLOCATE( svf(ndsvf,nsvfl) ) 6778 IF ( .NOT. ALLOCATED( svfsurf ) ) ALLOCATE( svfsurf(idsvf,nsvfl) ) 6779 READ(fsvf) svf 6780 READ(fsvf) svfsurf 6781 ENDIF 6782 6783 IF ( plant_canopy .AND. ncsfl > 0 ) THEN 6784 IF ( .NOT. ALLOCATED( csf ) ) ALLOCATE( csf(ndcsf,ncsfl) ) 6785 IF ( .NOT. ALLOCATED( csfsurf ) ) ALLOCATE( csfsurf(idcsf,ncsfl) ) 6786 READ(fsvf) csf 6787 READ(fsvf) csfsurf 6763 6788 ENDIF 6764 6789 READ ( fsvf ) svf_code_field 6765 6790 6766 6791 IF ( TRIM(svf_code_field) /= TRIM(svf_code) ) THEN 6767 6768 6792 WRITE( message_string, * ) 'Wrong structure of binary svf file' 6793 CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 ) 6769 6794 ENDIF 6770 6795 ! 6771 !-- Close binary file6772 CALL close_file( fsvf )6796 !-- Close binary file 6797 CALL close_file( fsvf ) 6773 6798 6774 6799 ENDIF 6775 6800 #if defined( __parallel ) 6776 6801 CALL MPI_BARRIER( comm2d, ierr ) 6777 6802 #endif 6778 6803 ENDDO … … 6801 6826 WRITE ( fsvf ) rad_version 6802 6827 WRITE ( fsvf ) nsvfl, ncsfl, nsurfl 6803 WRITE ( fsvf ) skyvf 6804 WRITE ( fsvf ) skyvft 6805 WRITE ( fsvf ) svf 6806 WRITE ( fsvf ) svfsurf 6807 IF ( plant_canopy ) THEN 6828 IF ( nsurfl > 0 ) THEN 6829 WRITE ( fsvf ) skyvf 6830 WRITE ( fsvf ) skyvft 6831 ENDIF 6832 IF ( nsvfl > 0 ) THEN 6833 WRITE ( fsvf ) svf 6834 WRITE ( fsvf ) svfsurf 6835 ENDIF 6836 IF ( plant_canopy .AND. ncsfl > 0 ) THEN 6808 6837 WRITE ( fsvf ) csf 6809 6838 WRITE ( fsvf ) csfsurf -
palm/trunk/SOURCE/surface_mod.f90
r2970 r2977 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Implement changes from branch radiation (r2948-2971) with minor modifications, 29 ! plus some formatting. 30 ! (moh.hefny) 31 ! Added flag to check the existence of vertical urban/land surfaces, required 32 ! to activate RTM 33 ! 34 ! 2970 2018-04-13 15:09:23Z suehring 28 35 ! Remove un-necessary initialization of surface elements in old large-scale 29 36 ! forcing mode … … 450 457 INTEGER(iwp) :: ns_v_on_file(0:3) !< total number of vertical surfaces with the same facing, required for writing restart data 451 458 459 LOGICAL :: vertical_surfaces_exist = .FALSE. !< flag indicating that there are vertical urban/land surfaces 460 !< in the domain (required to activiate RTM) 461 452 462 453 463 SAVE … … 488 498 PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, & 489 499 ns_v_on_file, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & 490 surf_usm_h, surf_usm_v, surf_type 500 surf_usm_h, surf_usm_v, surf_type, vertical_surfaces_exist 491 501 ! 492 502 !-- Public subroutines and functions … … 615 625 SUBROUTINE init_surface_arrays 616 626 627 628 USE pegrid 629 630 617 631 IMPLICIT NONE 618 632 … … 628 642 INTEGER(iwp), DIMENSION(0:3) :: num_lsm_v !< number of vertically-aligned natural surfaces 629 643 INTEGER(iwp), DIMENSION(0:3) :: num_usm_v !< number of vertically-aligned urban surfaces 644 645 INTEGER(iwp) :: num_surf_v_l !< number of vertically-aligned local urban/land surfaces 646 INTEGER(iwp) :: num_surf_v !< number of vertically-aligned total urban/land surfaces 630 647 631 648 LOGICAL :: building !< flag indicating building grid point … … 920 937 nys, nyn, nxl, nxr ) 921 938 ENDDO 922 939 ! 940 !-- Set the flag for the existence of vertical urban/land surfaces 941 num_surf_v_l = 0 942 DO l = 0, 3 943 num_surf_v_l = num_surf_v_l + surf_usm_v(l)%ns + surf_lsm_v(l)%ns 944 ENDDO 945 946 #if defined( __parallel ) 947 CALL MPI_ALLREDUCE( num_surf_v_l, num_surf_v, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr) 948 #else 949 num_surf_v = num_surf_v_l 950 #endif 951 IF ( num_surf_v > 0 ) vertical_surfaces_exist = .TRUE. 952 923 953 924 954 END SUBROUTINE init_surface_arrays -
palm/trunk/SOURCE/time_integration.f90
r2941 r2977 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implement changes from branch radiation (r2948-2971) with minor modifications. 28 ! (moh.hefny) 29 ! Fixed bug in if statement 30 ! 31 ! 2941 2018-04-03 11:54:58Z kanani 27 32 ! Deduct spinup_time from RUN_CONTROL output of main 3d run 28 33 ! (use time_since_reference_point instead of simulated_time) … … 1034 1039 CALL cpu_log( log_point(50), 'radiation', 'stop' ) 1035 1040 1036 IF ( urban_surface .OR. land_surface .AND.&1041 IF ( ( urban_surface .OR. land_surface ) .AND. & 1037 1042 radiation_interactions ) THEN 1038 1043 CALL cpu_log( log_point(75), 'radiation_interaction', 'start' ) -
palm/trunk/SOURCE/urban_surface_mod.f90
r2963 r2977 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Implement changes from branch radiation (r2948-2971) with minor modifications. 31 ! (moh.hefny): 32 ! Extended exn for all model domain height to avoid the need to get nzut. 33 ! 34 ! 2963 2018-04-12 14:47:44Z suehring 30 35 ! Introduce index for vegetation/wall, pavement/green-wall and water/window 31 36 ! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. . … … 3303 3308 DO i = nxl, nxr 3304 3309 DO k = nzb_do, nzt_do 3305 ! print*, j,nys,nyn,i,nxl,nxr,k,nzb_do,nzt_do,local_pf(i,j,k),temp_pf(k,j,i)3306 3310 local_pf(i,j,k) = temp_pf(k,j,i) 3307 3311 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
