Changeset 3226 for palm/trunk/SOURCE


Ignore:
Timestamp:
Aug 31, 2018 12:27:09 PM (6 years ago)
Author:
suehring
Message:

Bugfixes in calculation of sky-view and canopy-sink factors

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3186 r3226  
    2828! -----------------
    2929! $Id$
     30! Bugfixes in calculation of sky-view factors and canopy-sink factors.
     31!
     32! 3186 2018-07-30 17:07:14Z suehring
    3033! Remove print statement
    3134!
     
    393396!> @todo Output of full column vertical profiles used in RRTMG
    394397!> @todo Output of other rrtm arrays (such as volume mixing ratios)
    395 !> @todo Adapt for use with topography
     398!> @todo Check for mis-used NINT() calls in raytrace_2d
    396399!> @todo Optimize radiation_tendency routines
    397400!>
     
    29993002          rrtm_tsfc = t_rad_urb
    30003003         
    3001           IF ( lw_radiation )  THEN
     3004          IF ( lw_radiation )  THEN       
     3005         
    30023006             CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
    30033007             rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
     
    55455549        REAL(wp)                                      :: horizon       !< computed horizon height (tangent of elevation)
    55465550        REAL(wp)                                      :: azen          !< zenith angle
     5551        REAL(wp), DIMENSION(2)                        :: yxdir         !< y,x *unit* vector of ray direction (in grid units)
    55475552        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zdirs         !< directions in z (tangent of elevation)
    55485553        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zbdry         !< zenith angle boundaries
     
    58235828                  !--sum of vffrac for all iaz equals 1, verified
    58245829               ENDIF
    5825                CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,      &
    5826                                     surfstart(myid) + isurflt, facearea(td),  &
    5827                                     vffrac, .TRUE., .FALSE., win_lad, horizon,&
    5828                                     ztransp) !FIXME unit vect in grid units + zdirs
     5830               yxdir = (/ COS(azmid), SIN(azmid) /)
     5831               CALL raytrace_2d(ta, yxdir, zdirs,     &
     5832                                   -999, -999._wp, vffrac, .FALSE., .FALSE., &
     5833                                   win_lad, horizon, ztransp)
    58295834
    58305835               azen = pi/2 - ATAN(horizon)
     
    59515956           DO iaz = 1, naz
    59525957              azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
    5953               CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,     &
    5954                                    -999, -999._wp, vffrac, .FALSE., .TRUE., &
    5955                                    win_lad, horizon, ztransp) !FIXME unit vect in grid units + zdirs
     5958              yxdir = (/ COS(azmid), SIN(azmid) /)
     5959              CALL raytrace_2d(ta, yxdir, zdirs,                               &
     5960                               surfstart(myid) + isurflt, facearea(td),        &
     5961                               vffrac, .TRUE., .TRUE., win_lad, horizon,       &
     5962                               ztransp)
    59565963
    59575964              !--Save direct solar transparency
     
    65736580            ntrack = ntrack + 1
    65746581            crmid = (lastdist + nextdist) * .5_wp
    6575             column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
     6582            column = INT(yxorigin(:) + yxdir(:) * crmid, iwp)  !NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
    65766583
    65776584!--         calculate index of the grid with global indices (column(1),column(2))
     
    66586665         !--
    66596666         maxboxes = (ntrack + MAX(origin(1) - nzub, nzpt - origin(1))) * SIZE(zdirs, 1)
    6660          IF ( ncsfl + maxboxes > ncsfla )  THEN
    6661 !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
    6662 !--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
    6663 !--                                                / log(grow_factor)), kind=wp))
    6664 !--         or use this code to simply always keep some extra space after growing
    6665             k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
    6666             CALL merge_and_grow_csf(k)
    6667          ENDIF
     6667!          IF ( ncsfl + maxboxes > ncsfla )  THEN
     6668! !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
     6669! !--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
     6670! !--                                                / log(grow_factor)), kind=wp))
     6671! !--         or use this code to simply always keep some extra space after growing
     6672!             k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
     6673!             CALL merge_and_grow_csf(k)
     6674!          ENDIF
    66686675
    66696676         !--Calculate transparencies and store new CSFs
     
    67106717                     IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    67116718                        curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     6719                        IF (ncsfl+1 > SIZE(acsf) )                             &
     6720                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
    67126721
    67136722                        ncsfl = ncsfl + 1
     
    67706779
    67716780                     IF ( create_csf )  THEN
     6781                        IF (ncsfl+1 > SIZE(acsf) )                             &
     6782                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
    67726783                        ncsfl = ncsfl + 1
    67736784                        acsf(ncsfl)%ip = ip
Note: See TracChangeset for help on using the changeset viewer.