Changeset 4227 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- Sep 10, 2019 6:04:34 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r4226 r4227 28 28 ! ----------------- 29 29 ! $Id$ 30 ! implement new palm_date_time_mod 31 ! 32 ! 4226 2019-09-10 17:03:24Z suehring 30 33 ! - Netcdf input routine for dimension length renamed 31 34 ! - Define time variable for external radiation input relative to time_utc_init … … 232 235 ONLY: ddx, ddy, dx, dy 233 236 234 USE date_and_time_mod, &235 ONLY: calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,&236 day_of_year, d_seconds_year, day_of_month, day_of_year_init, &237 init_date_and_time, month_of_year, time_utc_init, time_utc238 239 237 USE indices, & 240 238 ONLY: nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, & … … 275 273 real_1d_3d, & 276 274 vars_pids 275 276 USE palm_date_time_mod, & 277 ONLY: date_time_str_len, get_date_time, & 278 hours_per_day, seconds_per_hour 277 279 278 280 USE plant_canopy_model_mod, & … … 387 389 time_radiation = 0.0_wp !< time since last call of radiation code 388 390 391 INTEGER(iwp) :: day_of_year !< day of the current year 389 392 390 393 REAL(wp) :: cos_zenith !< cosine of solar zenith angle, also z-coordinate of solar unit vector 394 REAL(wp) :: d_hours_day !< 1 / hours-per-day 395 REAL(wp) :: d_seconds_hour !< 1 / seconds-per-hour 396 REAL(wp) :: second_of_day !< second of the current day 391 397 REAL(wp) :: sun_dir_lat !< y-coordinate of solar unit vector 392 398 REAL(wp) :: sun_dir_lon !< x-coordinate of solar unit vector … … 1520 1526 CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 ) 1521 1527 ENDIF 1528 ! 1529 !-- Precalculate some time constants 1530 d_hours_day = 1.0_wp / REAL( hours_per_day, KIND = wp ) 1531 d_seconds_hour = 1.0_wp / seconds_per_hour 1532 1522 1533 ! 1523 1534 !-- If required, initialize radiation interactions between surfaces … … 2280 2291 !-- Calculate initial values of current (cosine of) the zenith angle and 2281 2292 !-- whether the sun is up 2282 CALL calc_zenith2283 ! 2284 !-- readjust date and time to its initial value 2285 CALL init_date_and_time2293 CALL get_date_time( time_since_reference_point, & 2294 day_of_year=day_of_year, & 2295 second_of_day=second_of_day ) 2296 CALL calc_zenith( day_of_year, second_of_day ) 2286 2297 ! 2287 2298 !-- Calculate initial surface albedo for different surfaces … … 2640 2651 ENDIF 2641 2652 2642 IF ( time_rad_f%var1d(0) /= time_utc_init ) THEN 2653 CALL get_date_time( 0.0_wp, second_of_day=second_of_day ) 2654 2655 IF ( time_rad_f%var1d(0) /= second_of_day ) THEN 2643 2656 message_string = 'External radiation forcing: first point in ' // & 2644 'time is /= time_utc_init.'2657 'time is /= origin_date_time.' 2645 2658 CALL message( 'radiation_init', 'PA0313', 1, 2, 0, 6, 0 ) 2646 2659 ENDIF 2647 2660 2648 2661 IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) & 2649 - time_utc_init) THEN2662 - second_of_day ) THEN 2650 2663 message_string = 'External radiation forcing does not cover ' // & 2651 2664 'the entire simulation time.' … … 2765 2778 2766 2779 END SELECT 2767 !2768 !-- Readjust date and time to its initial value2769 CALL init_date_and_time2770 2780 2771 2781 ! … … 2821 2831 LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces 2822 2832 2823 REAL(wp) :: fac_dt !< interpolation factor 2833 REAL(wp) :: fac_dt !< interpolation factor 2834 REAL(wp) :: second_of_day_init !< second of the day at model start 2824 2835 2825 2836 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine … … 2827 2838 ! 2828 2839 !-- Calculate current zenith angle 2829 CALL calc_zenith 2840 CALL get_date_time( time_since_reference_point, & 2841 day_of_year=day_of_year, & 2842 second_of_day=second_of_day ) 2843 CALL calc_zenith( day_of_year, second_of_day ) 2830 2844 ! 2831 2845 !-- Interpolate external radiation on current timestep … … 2835 2849 fac_dt = 0 2836 2850 ELSE 2851 CALL get_date_time( 0.0_wp, second_of_day=second_of_day_init ) 2837 2852 t = 0 2838 2853 DO WHILE ( time_rad_f%var1d(t) <= & 2839 time_since_reference_point + time_utc_init )2854 time_since_reference_point + second_of_day_init ) 2840 2855 t = t + 1 2841 2856 ENDDO … … 2843 2858 tm = MAX( t-1, 0 ) 2844 2859 2845 fac_dt = ( time_since_reference_point + time_utc_init&2860 fac_dt = ( time_since_reference_point + second_of_day_init & 2846 2861 - time_rad_f%var1d(tm) + dt_3d ) & 2847 2862 / ( time_rad_f%var1d(t) - time_rad_f%var1d(tm) ) … … 3067 3082 SUBROUTINE radiation_clearsky 3068 3083 3069 3070 3084 IMPLICIT NONE 3071 3085 3072 3086 INTEGER(iwp) :: l !< running index for surface orientation 3073 3087 3074 3088 LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces 3075 3089 3076 3090 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 3077 3091 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain … … 3083 3097 ! 3084 3098 !-- Calculate current zenith angle 3085 CALL calc_zenith 3099 CALL get_date_time( time_since_reference_point, & 3100 day_of_year=day_of_year, & 3101 second_of_day=second_of_day ) 3102 CALL calc_zenith( day_of_year, second_of_day ) 3086 3103 3087 3104 ! … … 3632 3649 ONLY: nbgp 3633 3650 3651 USE palm_date_time_mod, & 3652 ONLY: hours_per_day 3653 3634 3654 USE particle_attributes, & 3635 3655 ONLY: grid_particles, number_of_particles, particles, prt_count … … 3642 3662 INTEGER(iwp) :: k_topo !< topography top index 3643 3663 3664 REAL(wp) :: d_hours_day !< 1 / hours-per-day 3644 3665 REAL(wp) :: nc_rad, & !< number concentration of cloud droplets 3645 3666 s_r2, & !< weighted sum over all droplets with r^2 … … 3662 3683 3663 3684 ! 3685 !-- Pre-calculate parameters 3686 d_hours_day = 1.0_wp / REAL( hours_per_day, KIND=wp ) 3687 3688 ! 3664 3689 !-- Calculate current (cosine of) zenith angle and whether the sun is up 3665 CALL calc_zenith 3690 CALL get_date_time( time_since_reference_point, & 3691 day_of_year=day_of_year, & 3692 second_of_day=second_of_day ) 3693 CALL calc_zenith( day_of_year, second_of_day ) 3666 3694 zenith(0) = cos_zenith 3667 3695 ! … … 4548 4576 !> Calculate the cosine of the zenith angle (variable is called zenith) 4549 4577 !------------------------------------------------------------------------------! 4550 SUBROUTINE calc_zenith 4578 SUBROUTINE calc_zenith( day_of_year, second_of_day ) 4579 4580 USE palm_date_time_mod, & 4581 ONLY: seconds_per_day 4551 4582 4552 4583 IMPLICIT NONE 4553 4584 4554 REAL(wp) :: declination, & !< solar declination angle 4555 hour_angle !< solar hour angle 4556 ! 4557 !-- Calculate current day and time based on the initial values and simulation 4558 !-- time 4559 CALL calc_date_and_time 4560 4561 ! 4562 !-- Calculate solar declination and hour angle 4585 INTEGER(iwp), INTENT(IN) :: day_of_year !< day of the year 4586 4587 REAL(wp) :: declination !< solar declination angle 4588 REAL(wp) :: hour_angle !< solar hour angle 4589 REAL(wp), INTENT(IN) :: second_of_day !< current time of the day in UTC 4590 4591 ! 4592 !-- Calculate solar declination and hour angle 4563 4593 declination = ASIN( decl_1 * SIN(decl_2 * REAL(day_of_year, KIND=wp) - decl_3) ) 4564 hour_angle = 2.0_wp * pi * ( time_utc / 86400.0_wp) + lon - pi4594 hour_angle = 2.0_wp * pi * ( second_of_day / seconds_per_day ) + lon - pi 4565 4595 4566 4596 ! … … 5417 5447 5418 5448 sun_direction = .TRUE. 5419 CALL calc_zenith !< required also for diffusion radiation 5420 5449 CALL get_date_time( time_since_reference_point, & 5450 day_of_year=day_of_year, & 5451 second_of_day=second_of_day ) 5452 CALL calc_zenith( day_of_year, second_of_day ) !< required also for diffusion radiation 5453 5454 ! 5421 5455 !-- prepare rotated normal vectors and irradiance factor 5422 5456 vnorm(1,:) = kdir(:) … … 6240 6274 !------------------------------------------------------------------------------! 6241 6275 SUBROUTINE calc_diffusion_radiation 6242 6243 INTEGER(iwp) :: i !< grid index x-direction 6244 INTEGER(iwp) :: j !< grid index y-direction 6245 6246 REAL(wp) :: year_angle !< angle 6247 REAL(wp) :: etr !< extraterestrial radiation 6248 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 6249 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 6250 REAL(wp) :: clearnessIndex !< clearness index 6251 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 6252 6253 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 6254 ! 6276 6277 USE palm_date_time_mod, & 6278 ONLY: seconds_per_day 6279 6280 INTEGER(iwp) :: i !< grid index x-direction 6281 INTEGER(iwp) :: j !< grid index y-direction 6282 INTEGER(iwp) :: days_per_year !< days in the current year 6283 6284 REAL(wp) :: clearnessIndex !< clearness index 6285 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 6286 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 6287 REAL(wp) :: etr !< extraterestrial radiation 6288 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 6289 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 6290 REAL(wp) :: second_of_year !< current second of the year 6291 REAL(wp) :: year_angle !< angle 6292 6293 ! 6255 6294 !-- Calculate current day and time based on the initial values and simulation time 6256 year_angle = ( (day_of_year_init * 86400) + time_utc_init & 6257 + time_since_reference_point ) * d_seconds_year & 6258 * 2.0_wp * pi 6259 6295 CALL get_date_time( time_since_reference_point, & 6296 second_of_year = second_of_year, & 6297 days_per_year = days_per_year ) 6298 year_angle = second_of_year / ( REAL( days_per_year, KIND=wp ) * seconds_per_day ) & 6299 * 2.0_wp * pi 6300 6260 6301 etr = solar_constant * (1.00011_wp + & 6261 6302 0.034221_wp * cos(year_angle) + & … … 6263 6304 0.000719_wp * cos(2.0_wp * year_angle) + & 6264 6305 0.000077_wp * sin(2.0_wp * year_angle)) 6265 !6306 6266 6307 !-- 6267 6308 !-- Under a very low angle, we keep extraterestrial radiation at … … 8564 8605 IMPLICIT NONE 8565 8606 8566 INTEGER(iwp) :: it, i, j 8567 INTEGER(iwp) :: day_of_month_prev,month_of_year_prev 8568 REAL(wp) :: tsrp_prev 8569 REAL(wp) :: simulated_time_prev 8570 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsidir_tmp !< dsidir_tmp[:,i] = unit vector of i-th 8571 !< appreant solar direction 8607 INTEGER(iwp) :: it, i, j !< loop indices 8608 8609 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsidir_tmp !< dsidir_tmp[:,i] = unit vector of i-th 8610 !< appreant solar direction 8572 8611 8573 8612 ALLOCATE ( dsidir_rev(0:raytrace_discrete_elevs/2-1, & … … 8577 8616 raytrace_discrete_elevs/2*raytrace_discrete_azims) ) 8578 8617 ndsidir = 0 8579 8580 !8581 !-- We will artificialy update time_since_reference_point and return to8582 !-- true value later8583 tsrp_prev = time_since_reference_point8584 simulated_time_prev = simulated_time8585 day_of_month_prev = day_of_month8586 month_of_year_prev = month_of_year8587 8618 sun_direction = .TRUE. 8588 8589 ! 8590 !-- initialize the simulated_time 8591 simulated_time = 0._wp 8619 8592 8620 ! 8593 8621 !-- Process spinup time if configured 8594 8622 IF ( spinup_time > 0._wp ) THEN 8595 8623 DO it = 0, CEILING(spinup_time / dt_spinup) 8596 time_since_reference_point = -spinup_time + REAL(it, wp) * dt_spinup 8597 simulated_time = simulated_time + dt_spinup 8598 CALL simulate_pos 8624 CALL simulate_pos( it * dt_spinup - spinup_time ) 8599 8625 ENDDO 8600 8626 ENDIF … … 8602 8628 !-- Process simulation time 8603 8629 DO it = 0, CEILING(( end_time - spinup_time ) / dt_radiation) 8604 time_since_reference_point = REAL(it, wp) * dt_radiation 8605 simulated_time = simulated_time + dt_radiation 8606 CALL simulate_pos 8630 CALL simulate_pos( it * dt_radiation ) 8607 8631 ENDDO 8608 8632 ! 8609 !-- Return date and time to its original values8610 time_since_reference_point = tsrp_prev8611 simulated_time = simulated_time_prev8612 day_of_month = day_of_month_prev8613 month_of_year = month_of_year_prev8614 CALL init_date_and_time8615 8616 8633 !-- Allocate global vars which depend on ndsidir 8617 8634 ALLOCATE ( dsidir ( 3, ndsidir ) ) … … 8634 8651 !> Simuates a single position 8635 8652 !------------------------------------------------------------------------! 8636 SUBROUTINE simulate_pos 8637 IMPLICIT NONE 8653 SUBROUTINE simulate_pos( time_since_reference_local ) 8654 8655 REAL(wp), INTENT(IN) :: time_since_reference_local !< local time since reference 8638 8656 ! 8639 8657 !-- Update apparent solar position based on modified t_s_r_p 8640 CALL calc_zenith 8658 CALL get_date_time( time_since_reference_local, & 8659 day_of_year=day_of_year, & 8660 second_of_day=second_of_day ) 8661 CALL calc_zenith( day_of_year, second_of_day ) 8641 8662 IF ( cos_zenith > 0 ) THEN 8642 8663 !--
Note: See TracChangeset
for help on using the changeset viewer.
