Changeset 3226 for palm/trunk/SOURCE

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:

• palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (10 diffs)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Bugfixes in calculation of sky-view factors and canopy-sink factors.
+! 
+! 3186 2018-07-30 17:07:14Z suehring
 ! Remove print statement
 ! 
@@ -393 +396 @@
 !> @todo Output of full column vertical profiles used in RRTMG
 !> @todo Output of other rrtm arrays (such as volume mixing ratios)
-!> @todo Adapt for use with topography
+!> @todo Check for mis-used NINT() calls in raytrace_2d
 !> @todo Optimize radiation_tendency routines
 !>
@@ -2999 +3002 @@
           rrtm_tsfc = t_rad_urb
           
-          IF ( lw_radiation )  THEN
+          IF ( lw_radiation )  THEN        
+          
              CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
              rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
@@ -5545 +5549 @@
         REAL(wp)                                      :: horizon       !< computed horizon height (tangent of elevation)
         REAL(wp)                                      :: azen          !< zenith angle
+        REAL(wp), DIMENSION(2)                        :: yxdir         !< y,x *unit* vector of ray direction (in grid units)
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: zdirs         !< directions in z (tangent of elevation)
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: zbdry         !< zenith angle boundaries
@@ -5823 +5828 @@
                   !--sum of vffrac for all iaz equals 1, verified
                ENDIF
-               CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,      &
-                                    surfstart(myid) + isurflt, facearea(td),  &
-                                    vffrac, .TRUE., .FALSE., win_lad, horizon,&
-                                    ztransp) !FIXME unit vect in grid units + zdirs
+               yxdir = (/ COS(azmid), SIN(azmid) /)
+               CALL raytrace_2d(ta, yxdir, zdirs,     &
+                                   -999, -999._wp, vffrac, .FALSE., .FALSE., &
+                                   win_lad, horizon, ztransp)
 
                azen = pi/2 - ATAN(horizon)
@@ -5951 +5956 @@
            DO iaz = 1, naz
               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
-              CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,     &
-                                   -999, -999._wp, vffrac, .FALSE., .TRUE., &
-                                   win_lad, horizon, ztransp) !FIXME unit vect in grid units + zdirs
+              yxdir = (/ COS(azmid), SIN(azmid) /)
+              CALL raytrace_2d(ta, yxdir, zdirs,                               &
+                               surfstart(myid) + isurflt, facearea(td),        &
+                               vffrac, .TRUE., .TRUE., win_lad, horizon,       &
+                               ztransp)
 
               !--Save direct solar transparency
@@ -6573 +6580 @@
             ntrack = ntrack + 1
             crmid = (lastdist + nextdist) * .5_wp
-            column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
+            column = INT(yxorigin(:) + yxdir(:) * crmid, iwp)  !NINT(yxorigin(:) + yxdir(:) * crmid, iwp) 
 
 !--         calculate index of the grid with global indices (column(1),column(2))
@@ -6658 +6665 @@
          !--
          maxboxes = (ntrack + MAX(origin(1) - nzub, nzpt - origin(1))) * SIZE(zdirs, 1)
-         IF ( ncsfl + maxboxes > ncsfla )  THEN
-!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
-!--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
-!--                                                / log(grow_factor)), kind=wp))
-!--         or use this code to simply always keep some extra space after growing
-            k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
-            CALL merge_and_grow_csf(k)
-         ENDIF
+!          IF ( ncsfl + maxboxes > ncsfla )  THEN
+! !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
+! !--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
+! !--                                                / log(grow_factor)), kind=wp))
+! !--         or use this code to simply always keep some extra space after growing
+!             k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
+!             CALL merge_and_grow_csf(k)
+!          ENDIF
 
          !--Calculate transparencies and store new CSFs
@@ -6710 +6717 @@
                      IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
                         curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
+                        IF (ncsfl+1 > SIZE(acsf) )                             &
+                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
 
                         ncsfl = ncsfl + 1
@@ -6770 +6779 @@
 
                      IF ( create_csf )  THEN
+                        IF (ncsfl+1 > SIZE(acsf) )                             &
+                           CALL merge_and_grow_csf(CEILING(REAL(ncsfl, kind=wp) * grow_factor))
                         ncsfl = ncsfl + 1
                         acsf(ncsfl)%ip = ip