Changeset 3272 for palm/trunk
- Timestamp:
- Sep 24, 2018 10:16:32 AM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90
r1585 r3272 78 78 subroutine rrtmg_sw & 79 79 (ncol ,nlay ,icld ,iaer , & 80 play ,plev ,tlay ,tlev ,tsfc , & 81 h2ovmr ,o3vmr ,co2vmr ,ch4vmr ,n2ovmr ,o2vmr, & 82 asdir ,asdif ,aldir ,aldif , & 83 coszen ,adjes ,dyofyr ,scon , & 84 inflgsw ,iceflgsw,liqflgsw,cldfr , & 85 taucld ,ssacld ,asmcld ,fsfcld , & 86 cicewp ,cliqwp ,reice ,reliq , & 87 tauaer ,ssaaer ,asmaer ,ecaer , & 88 swuflx ,swdflx ,swhr ,swuflxc ,swdflxc ,swhrc) 80 play ,plev ,tlay ,tlev , & 81 tsfc ,h2ovmr ,o3vmr ,co2vmr , & 82 ch4vmr ,n2ovmr ,o2vmr ,asdir , & 83 asdif ,aldir ,aldif ,coszen , & 84 adjes ,dyofyr ,scon ,inflgsw , & 85 iceflgsw,liqflgsw,cldfr ,taucld , & 86 ssacld ,asmcld ,fsfcld ,cicewp , & 87 cliqwp ,reice ,reliq ,tauaer , & 88 ssaaer ,asmaer ,ecaer ,swuflx , & 89 swdflx ,swhr ,swuflxc ,swdflxc , & 90 swhrc ,dirdflux,difdflux) 89 91 90 92 ! ------- Description ------- … … 169 171 ! delta scaling based on setting of idelm flag. 170 172 ! Dec 2008: M. J. Iacono, AER, Inc. 173 !-- Added the output direct and diffuse fluxes to the argument of rrtmg_sw 174 ! Aug 2018: M. Salim, Humboldt Uni, Berlin, Germany 171 175 172 176 ! --------- Modules --------- … … 417 421 real(kind=rb) :: swnflx(nlay+2) ! Total sky shortwave net flux (W/m2) 418 422 real(kind=rb) :: swnflxc(nlay+2) ! Clear sky shortwave net flux (W/m2) 419 real(kind=rb) :: dirdflux(nlay+2) ! Direct downward shortwave surface flux420 real(kind=rb) :: difdflux(nlay+2) ! Diffuse downward shortwave surface flux423 real(kind=rb), intent(out) :: dirdflux(nlay+2) ! Direct downward shortwave surface flux 424 real(kind=rb), intent(out) :: difdflux(nlay+2) ! Diffuse downward shortwave surface flux 421 425 real(kind=rb) :: uvdflx(nlay+2) ! Total sky downward shortwave flux, UV/vis 422 426 real(kind=rb) :: nidflx(nlay+2) ! Total sky downward shortwave flux, near-IR -
palm/trunk/SOURCE/radiation_model_mod.f90
r3261 r3272 28 28 ! ----------------- 29 29 ! $Id$ 30 ! - split direct and diffusion shortwave radiation using RRTMG rather than using 31 ! calc_diffusion_radiation, in case of RRTMG 32 ! - removed the namelist variable split_diffusion_radiation. Now splitting depends 33 ! on the choise of radiation radiation scheme 34 ! - removed calculating the rdiation flux for surfaces at the radiation scheme 35 ! in case of using RTM since it will be calculated anyway in the radiation 36 ! interaction routine. 37 ! - set SW radiation flux for surfaces to zero at night in case of no RTM is used 38 ! - precalculate the unit vector yxdir of ray direction to avoid the temporarly 39 ! array allocation during the subroutine call 40 ! - fixed a bug in calculating the max number of boxes ray can cross in the domain 41 ! 42 ! 3264 2018-09-20 13:54:11Z moh.hefny 30 43 ! Bugfix in raytrace_2d calls 31 44 ! … … 722 735 rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2) 723 736 rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d) 724 rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 725 737 rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 738 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m²) 739 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m²) 726 740 727 741 REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation … … 822 836 823 837 !-- configuration parameters (they can be setup in PALM config) 824 LOGICAL :: split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation825 !< (.F. in case the radiation model already does it)826 838 LOGICAL :: rma_lad_raytrace = .FALSE. !< use MPI RMA to access LAD for raytracing (instead of global array) 827 839 LOGICAL :: mrt_factors = .FALSE. !< whether to generate MRT factor files during init … … 1080 1092 skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,& 1081 1093 zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, & 1082 split_diffusion_radiation, &1083 1094 nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf, & 1084 1095 svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir, & … … 1153 1164 SELECT CASE ( TRIM( var ) ) 1154 1165 1155 CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )1166 CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' ) 1156 1167 IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN 1157 1168 message_string = '"output of "' // TRIM( var ) // '" requi' // & … … 1189 1200 IF ( TRIM( var ) == 'rad_sw_in*' ) unit = 'W/m2' 1190 1201 IF ( TRIM( var ) == 'rad_sw_out*' ) unit = 'W/m2' 1202 IF ( TRIM( var ) == 'rad_sw_in' ) unit = 'W/m2' 1191 1203 IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = '' 1192 1204 IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = '' … … 2769 2781 radiation_scheme, skip_time_do_radiation, & 2770 2782 sw_radiation, unscheduled_radiation_calls, & 2771 split_diffusion_radiation, &2772 2783 max_raytracing_dist, min_irrf_value, & 2773 2784 nrefsteps, mrt_factors, rma_lad_raytrace, & … … 2782 2793 radiation_scheme, skip_time_do_radiation, & 2783 2794 sw_radiation, unscheduled_radiation_calls, & 2784 split_diffusion_radiation, &2785 2795 max_raytracing_dist, min_irrf_value, & 2786 2796 nrefsteps, mrt_factors, rma_lad_raytrace, & … … 3054 3064 rad_lw_out(k,:,:) = rrtm_lwuflx(0,k) 3055 3065 ENDDO 3056 3066 rad_lw_in_diff(:,:) = rad_lw_in(0,:,:) 3057 3067 ! 3058 3068 !-- Save heating rates (convert from K/d to K/h). … … 3070 3080 ENDDO 3071 3081 3072 !3073 !-- Save surface radiative fluxes and change in LW heating rate3074 !-- onto respective surface elements3075 !-- Horizontal surfaces3076 IF ( surf_lsm_h%ns > 0 ) THEN3077 surf_lsm_h%rad_lw_in = rrtm_lwdflx(0,nzb)3078 surf_lsm_h%rad_lw_out = rrtm_lwuflx(0,nzb)3079 surf_lsm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)3080 ENDIF3081 IF ( surf_usm_h%ns > 0 ) THEN3082 surf_usm_h%rad_lw_in = rrtm_lwdflx(0,nzb)3083 surf_usm_h%rad_lw_out = rrtm_lwuflx(0,nzb)3084 surf_usm_h%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)3085 ENDIF3086 !3087 !-- Vertical surfaces.3088 DO l = 0, 33089 IF ( surf_lsm_v(l)%ns > 0 ) THEN3090 surf_lsm_v(l)%rad_lw_in = rrtm_lwdflx(0,nzb)3091 surf_lsm_v(l)%rad_lw_out = rrtm_lwuflx(0,nzb)3092 surf_lsm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)3093 ENDIF3094 IF ( surf_usm_v(l)%ns > 0 ) THEN3095 surf_usm_v(l)%rad_lw_in = rrtm_lwdflx(0,nzb)3096 surf_usm_v(l)%rad_lw_out = rrtm_lwuflx(0,nzb)3097 surf_usm_v(l)%rad_lw_out_change_0 = rrtm_lwuflx_dt(0,nzb)3098 ENDIF3099 ENDDO3100 3101 3082 ENDIF 3102 3083 3103 3084 IF ( sw_radiation .AND. sun_up ) THEN 3104 CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer,&3105 rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev,&3106 rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr,&3107 rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir,&3108 rrtm_asdif , rrtm_aldir , rrtm_aldif , zenith,&3109 0.0_wp , day_of_year , solar_constant, rrtm_inflgsw,&3110 rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld,&3111 rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp,&3112 rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer,&3113 rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer ,&3114 rrtm_sw uflx , rrtm_swdflx , rrtm_swhr ,&3115 rrtm_sw uflxc , rrtm_swdflxc , rrtm_swhrc)3085 CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,& 3086 rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,& 3087 rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,& 3088 rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir ,& 3089 rrtm_asdif , rrtm_aldir , rrtm_aldif , zenith ,& 3090 0.0_wp , day_of_year , solar_constant, rrtm_inflgsw ,& 3091 rrtm_iceflgsw , rrtm_liqflgsw , rrtm_cldfr , rrtm_sw_taucld ,& 3092 rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,& 3093 rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,& 3094 rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx ,& 3095 rrtm_swdflx , rrtm_swhr , rrtm_swuflxc , rrtm_swdflxc ,& 3096 rrtm_swhrc , rrtm_dirdflux , rrtm_difdflux ) 3116 3097 3117 3098 ! 3118 !-- Save fluxes 3099 !-- Save fluxes: 3100 !-- - whole domain 3119 3101 DO k = nzb, nzt+1 3120 3102 rad_sw_in(k,:,:) = rrtm_swdflx(0,k) 3121 3103 rad_sw_out(k,:,:) = rrtm_swuflx(0,k) 3122 3104 ENDDO 3105 !-- - direct and diffuse SW at urban-surface-layer (required by RTM) 3106 rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb) 3107 rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb) 3123 3108 3124 3109 ! … … 3128 3113 rad_sw_cs_hr(k,:,:) = rrtm_swhrc(0,k) * d_hours_day 3129 3114 ENDDO 3130 3131 !3132 !-- Save surface radiative fluxes onto respective surface elements3133 !-- Horizontal surfaces3134 IF ( surf_lsm_h%ns > 0 ) THEN3135 surf_lsm_h%rad_sw_in = rrtm_swdflx(0,nzb)3136 surf_lsm_h%rad_sw_out = rrtm_swuflx(0,nzb)3137 ENDIF3138 IF ( surf_usm_h%ns > 0 ) THEN3139 surf_usm_h%rad_sw_in = rrtm_swdflx(0,nzb)3140 surf_usm_h%rad_sw_out = rrtm_swuflx(0,nzb)3141 ENDIF3142 !3143 !-- Vertical surfaces. Fluxes are obtain at respective vertical3144 !-- level of the surface element3145 DO l = 0, 33146 IF ( surf_lsm_v(l)%ns > 0 ) THEN3147 surf_lsm_v(l)%rad_sw_in = rrtm_swdflx(0,nzb)3148 surf_lsm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)3149 ENDIF3150 IF ( surf_usm_v(l)%ns > 0 ) THEN3151 surf_usm_v(l)%rad_sw_in = rrtm_swdflx(0,nzb)3152 surf_usm_v(l)%rad_sw_out = rrtm_swuflx(0,nzb)3153 ENDIF3154 ENDDO3155 3115 ! 3156 3116 !-- Solar radiation is zero during night … … 3158 3118 rad_sw_in = 0.0_wp 3159 3119 rad_sw_out = 0.0_wp 3120 rad_sw_in_dir(:,:) = 0.0_wp 3121 rad_sw_in_diff(:,:) = 0.0_wp 3160 3122 ENDIF 3161 3123 ! 3162 3124 !-- RRTMG is called for each (j,i) grid point separately, starting at the 3163 !-- highest topography level 3125 !-- highest topography level. Here no RTM is used since average_radiation is false 3164 3126 ELSE 3165 3127 ! … … 3551 3513 rrtm_swuflxc(:,k_topo:nzt_rad+1), & 3552 3514 rrtm_swdflxc(:,k_topo:nzt_rad+1), & 3553 rrtm_swhrc(:,k_topo+1:nzt_rad+1) ) 3515 rrtm_swhrc(:,k_topo+1:nzt_rad+1), & 3516 rrtm_dirdflux(:,k_topo:nzt_rad+1), & 3517 rrtm_difdflux(:,k_topo:nzt_rad+1) ) 3554 3518 3555 3519 DEALLOCATE( rrtm_sw_taucld_dum ) … … 3609 3573 rad_sw_in = 0.0_wp 3610 3574 rad_sw_out = 0.0_wp 3575 !-- !!!!!!!! ATTENSION !!!!!!!!!!!!!!! 3576 !-- Surface radiative fluxes should be also set to zero here 3577 !-- Save surface radiative fluxes onto respective surface elements 3578 !-- Horizontal surfaces 3579 DO m = surf_lsm_h%start_index(j,i), & 3580 surf_lsm_h%end_index(j,i) 3581 surf_lsm_h%rad_sw_in(m) = 0.0_wp 3582 surf_lsm_h%rad_sw_out(m) = 0.0_wp 3583 ENDDO 3584 DO m = surf_usm_h%start_index(j,i), & 3585 surf_usm_h%end_index(j,i) 3586 surf_usm_h%rad_sw_in(m) = 0.0_wp 3587 surf_usm_h%rad_sw_out(m) = 0.0_wp 3588 ENDDO 3589 ! 3590 !-- Vertical surfaces. Fluxes are obtain at respective vertical 3591 !-- level of the surface element 3592 DO l = 0, 3 3593 DO m = surf_lsm_v(l)%start_index(j,i), & 3594 surf_lsm_v(l)%end_index(j,i) 3595 k = surf_lsm_v(l)%k(m) 3596 surf_lsm_v(l)%rad_sw_in(m) = 0.0_wp 3597 surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp 3598 ENDDO 3599 DO m = surf_usm_v(l)%start_index(j,i), & 3600 surf_usm_v(l)%end_index(j,i) 3601 k = surf_usm_v(l)%k(m) 3602 surf_usm_v(l)%rad_sw_in(m) = 0.0_wp 3603 surf_usm_v(l)%rad_sw_out(m) = 0.0_wp 3604 ENDDO 3605 ENDDO 3611 3606 ENDIF 3612 3607 … … 3617 3612 ! 3618 3613 !-- Finally, calculate surface net radiation for surface elements. 3619 !-- First, for horizontal surfaces 3620 DO m = 1, surf_lsm_h%ns 3621 surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m) & 3622 - surf_lsm_h%rad_sw_out(m) & 3623 + surf_lsm_h%rad_lw_in(m) & 3624 - surf_lsm_h%rad_lw_out(m) 3625 ENDDO 3626 DO m = 1, surf_usm_h%ns 3627 surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m) & 3628 - surf_usm_h%rad_sw_out(m) & 3629 + surf_usm_h%rad_lw_in(m) & 3630 - surf_usm_h%rad_lw_out(m) 3631 ENDDO 3632 ! 3633 !-- Vertical surfaces. 3634 !-- Todo: weight with azimuth and zenith angle according to their orientation! 3635 DO l = 0, 3 3636 DO m = 1, surf_lsm_v(l)%ns 3637 surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m) & 3638 - surf_lsm_v(l)%rad_sw_out(m) & 3639 + surf_lsm_v(l)%rad_lw_in(m) & 3640 - surf_lsm_v(l)%rad_lw_out(m) 3614 IF ( .NOT. radiation_interactions ) THEN 3615 !-- First, for horizontal surfaces 3616 DO m = 1, surf_lsm_h%ns 3617 surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m) & 3618 - surf_lsm_h%rad_sw_out(m) & 3619 + surf_lsm_h%rad_lw_in(m) & 3620 - surf_lsm_h%rad_lw_out(m) 3641 3621 ENDDO 3642 DO m = 1, surf_usm_ v(l)%ns3643 surf_usm_ v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m)&3644 - surf_usm_v(l)%rad_sw_out(m)&3645 + surf_usm_v(l)%rad_lw_in(m)&3646 - surf_usm_v(l)%rad_lw_out(m)3622 DO m = 1, surf_usm_h%ns 3623 surf_usm_h%rad_net(m) = surf_usm_h%rad_sw_in(m) & 3624 - surf_usm_h%rad_sw_out(m) & 3625 + surf_usm_h%rad_lw_in(m) & 3626 - surf_usm_h%rad_lw_out(m) 3647 3627 ENDDO 3648 ENDDO 3628 ! 3629 !-- Vertical surfaces. 3630 !-- Todo: weight with azimuth and zenith angle according to their orientation! 3631 DO l = 0, 3 3632 DO m = 1, surf_lsm_v(l)%ns 3633 surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m) & 3634 - surf_lsm_v(l)%rad_sw_out(m) & 3635 + surf_lsm_v(l)%rad_lw_in(m) & 3636 - surf_lsm_v(l)%rad_lw_out(m) 3637 ENDDO 3638 DO m = 1, surf_usm_v(l)%ns 3639 surf_usm_v(l)%rad_net(m) = surf_usm_v(l)%rad_sw_in(m) & 3640 - surf_usm_v(l)%rad_sw_out(m) & 3641 + surf_usm_v(l)%rad_lw_in(m) & 3642 - surf_usm_v(l)%rad_lw_out(m) 3643 ENDDO 3644 ENDDO 3645 ENDIF 3649 3646 3650 3647 … … 3905 3902 DEALLOCATE ( rrtm_swhr ) 3906 3903 DEALLOCATE ( rrtm_swhrc ) 3904 DEALLOCATE ( rrtm_dirdflux ) 3905 DEALLOCATE ( rrtm_difdflux ) 3907 3906 3908 3907 ENDIF … … 4070 4069 ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) ) 4071 4070 ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) ) 4071 ALLOCATE ( rrtm_dirdflux(0:0,nzb:nzt_rad+1) ) 4072 ALLOCATE ( rrtm_difdflux(0:0,nzb:nzt_rad+1) ) 4072 4073 4073 4074 rrtm_swdflx = 0.0_wp … … 4077 4078 rrtm_swdflxc = 0.0_wp 4078 4079 rrtm_swhrc = 0.0_wp 4080 rrtm_dirdflux = 0.0_wp 4081 rrtm_difdflux = 0.0_wp 4079 4082 4080 4083 ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) ) … … 4482 4485 4483 4486 IMPLICIT NONE 4484 INTEGER(iwp) :: i, j, k, kk, d, refstep, m, mm, l, ll 4487 4488 INTEGER(iwp) :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll 4485 4489 INTEGER(iwp) :: isurf, isurfsrc, isvf, icsf, ipcgb 4486 4490 INTEGER(iwp) :: isd !< solar direction number … … 4519 4523 REAL(wp) :: area_hor !< total horizontal area of domain in all processor 4520 4524 4521 4522 4525 #if ! defined( __nopointer ) 4523 4526 IF ( plant_canopy ) THEN … … 4545 4548 4546 4549 IF ( zenith(0) > 0 ) THEN 4547 !-- 4548 !-- 4549 !-- 4550 !-- now we will "squash" the sunorig vector by grid box size in 4551 !-- each dimension, so that this new direction vector will allow us 4552 !-- to traverse the ray path within grid coordinates directly 4550 4553 sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /) 4551 !-- 4554 !-- sunorig_grid = sunorig_grid / norm2(sunorig_grid) 4552 4555 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) 4553 4556 4554 4557 IF ( npcbl > 0 ) THEN 4555 !-- 4558 !-- precompute effective box depth with prototype Leaf Area Density 4556 4559 pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1 4557 4560 CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift), & … … 4564 4567 ENDIF 4565 4568 4566 !-- split diffusion and direct part of the solar downward radiation 4567 !-- comming from radiation model and store it in 2D arrays 4568 !-- rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff 4569 IF ( split_diffusion_radiation ) THEN 4570 CALL calc_diffusion_radiation 4571 ELSE 4572 rad_sw_in_diff = 0.0_wp 4573 rad_sw_in_dir(:,:) = rad_sw_in(0,:,:) 4574 rad_lw_in_diff(:,:) = rad_lw_in(0,:,:) 4575 ENDIF 4569 !-- if radiation scheme is not RRTMG, split diffusion and direct part of the solar downward radiation 4570 !-- comming from radiation model and store it in 2D arrays 4571 IF ( radiation_scheme /= 'rrtmg' ) CALL calc_diffusion_radiation 4576 4572 4577 4573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 6696 6692 !-- 6697 6693 maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1) 6698 !IF ( ncsfl + maxboxes > ncsfla ) THEN6699 ! !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)6700 ! !-- k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &6701 ! !-- / log(grow_factor)), kind=wp))6702 ! !-- or use this code to simply always keep some extra space after growing6703 !k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)6704 !CALL merge_and_grow_csf(k)6705 !ENDIF6694 IF ( ncsfl + maxboxes > ncsfla ) THEN 6695 !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) 6696 !-- k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) & 6697 !-- / log(grow_factor)), kind=wp)) 6698 !-- or use this code to simply always keep some extra space after growing 6699 k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor) 6700 CALL merge_and_grow_csf(k) 6701 ENDIF 6706 6702 6707 6703 !--Calculate transparencies and store new CSFs … … 6748 6744 IF ( rt2_track_lad(iz, i) > 0._wp ) THEN 6749 6745 curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1))) 6750 IF (ncsfl+1 > SIZE(acsf) ) &6751 CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))6752 6753 6746 ncsfl = ncsfl + 1 6754 6747 acsf(ncsfl)%ip = ip … … 6810 6803 6811 6804 IF ( create_csf ) THEN 6812 IF (ncsfl+1 > SIZE(acsf) ) &6813 CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))6814 6805 ncsfl = ncsfl + 1 6815 6806 acsf(ncsfl)%ip = ip
Note: See TracChangeset
for help on using the changeset viewer.