Changeset 2209 for palm/trunk/SOURCE/urban_surface_mod.f90

Timestamp:

Apr 19, 2017 9:34:46 AM (7 years ago)

Author:

kanani

Message:

small bugfix, formatting and new PCM output

File:

• palm/trunk/SOURCE/urban_surface_mod.f90 (modified) (9 diffs)

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -21 +21 @@
 ! Current revisions:
 ! ------------------
-! 
+! cpp switch __mpi3 removed,
+! minor formatting,
+! small bugfix for division by zero (Krc)
 ! 
 ! Former revisions:
@@ -1489 +1491 @@
         CALL location_message( '    calculation of SVF and CSF', .TRUE. )
 
-#if defined( __mpi3 )
+
 !--     precalculate face areas for different face directions using normal vector
         DO d = 0, 9
@@ -1964 +1966 @@
         CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
        
-#endif
+
     END SUBROUTINE usm_calc_svf
 
@@ -2487 +2489 @@
 !------------------------------------------------------------------------------!
     PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)
-        IMPLICIT NONE
-        REAL(wp), DIMENSION(3), INTENT(in)      :: origin    !< ray origin
-        REAL(wp), DIMENSION(3), INTENT(in)      :: uvect     !< ray unit vector
-        INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross  !< found boundary crossing (d, z, y, x)
-        REAL(wp), DIMENSION(3)                  :: crossdist !< crossing distance
-        INTEGER(iwp), DIMENSION(3)              :: bdyd      !< boundary direction
-        REAL(wp)                                :: bdydim    !< 
-        REAL(wp)                                :: dist      !< 
-        INTEGER(iwp)                            :: seldim    !< found fist crossing index
-        INTEGER(iwp)                            :: d         !< 
-
-        bdydim = nzut + .5_wp !< top boundary
-        bdyd(1) = isky
-        crossdist(1) = (bdydim - origin(1)) / uvect(1)
-
-        IF ( uvect(2) >= 0._wp )  THEN
-            bdydim = ny + .5_wp !< north global boundary
-            bdyd(2) = inorthb
-        ELSE
-            bdydim = -.5_wp !< south global boundary
-            bdyd(2) = isouthb
-        ENDIF
-        crossdist(2) = (bdydim - origin(2)) / uvect(2)
-
-        IF ( uvect(3) >= 0._wp )  THEN
-            bdydim = nx + .5_wp !< east global boundary
-            bdyd(3) = ieastb
-        ELSE
-            bdydim = -.5_wp !< west global boundary
-            bdyd(3) = iwestb
-        ENDIF
-        crossdist(3) = (bdydim - origin(3)) / uvect(3)
-
-        seldim = minloc(crossdist, 1)
-        dist = crossdist(seldim)
-        d = bdyd(seldim)
-
-        bdycross(1) = d
-        bdycross(2:4) = NINT( origin(:) + uvect(:)*dist &
-                        + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
+    
+       IMPLICIT NONE
+       
+       INTEGER(iwp) ::  d       !<
+       INTEGER(iwp) ::  seldim  !< found fist crossing index
+
+       INTEGER(iwp), DIMENSION(3)              ::  bdyd      !< boundary direction       
+       INTEGER(iwp), DIMENSION(4), INTENT(out) ::  bdycross  !< found boundary crossing (d, z, y, x)
+       
+       REAL(wp)                                ::  bdydim  !< 
+       REAL(wp)                                ::  dist    !<
+       
+       REAL(wp), DIMENSION(3)             ::  crossdist  !< crossing distance
+       REAL(wp), DIMENSION(3), INTENT(in) ::  origin     !< ray origin
+       REAL(wp), DIMENSION(3), INTENT(in) ::  uvect      !< ray unit vector
+ 
+
+       bdydim       = nzut + .5_wp  !< top boundary
+       bdyd(1)      = isky
+       crossdist(1) = ( bdydim - origin(1) ) / uvect(1)  !< subroutine called only when uvect(1)>0
+
+       IF ( uvect(2) == 0._wp )  THEN
+          crossdist(2) = huge(1._wp)
+       ELSE
+          IF ( uvect(2) >= 0._wp )  THEN
+             bdydim  = ny + .5_wp  !< north global boundary
+             bdyd(2) = inorthb
+          ELSE
+             bdydim  = -.5_wp  !< south global boundary
+             bdyd(2) = isouthb
+          ENDIF
+          crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
+       ENDIF
+
+       IF ( uvect(3) == 0._wp )  THEN
+          crossdist(3) = huge(1._wp)
+       ELSE
+          IF ( uvect(3) >= 0._wp )  THEN
+             bdydim  = nx + .5_wp  !< east global boundary
+             bdyd(3) = ieastb
+          ELSE
+             bdydim  = -.5_wp  !< west global boundary
+             bdyd(3) = iwestb
+          ENDIF
+          crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
+       ENDIF
+
+       seldim = minloc(crossdist, 1)
+       dist   = crossdist(seldim)
+       d      = bdyd(seldim)
+
+       bdycross(1)   = d
+       bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
+                                       + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
+                        
     END SUBROUTINE
 
@@ -2864 +2881 @@
        urban_surface = .TRUE.
        
-!
-!--    Check whether pre-processor (cpp) option "__mpi3" is set. It is required
-!--    for the full functionality of the USM. "__mpi3" directive is implemented,
-!--    because some compilers cannot handle MPI-3 operations, hence, these parts
-!--    of code shall only be compiled if explicitly enabled.
-#if ! defined ( __mpi3 )
-          message_string = 'urban surface model requires compilation of ' //   & 
-                           'PALM with pre-processor directive -D__mpi3'
-          CALL message( 'usm_parin', 'PA0503', 1, 2, 0, 6, 0 )
-#endif
-
 
  10    CONTINUE
@@ -3220 +3226 @@
         REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
 
-#if defined( __mpi3 )
+
 !--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
 !--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
@@ -3363 +3369 @@
         
         visible = .TRUE.
-        
-#else
-        visible      = .FALSE.                          !Set variables to avoid compiler warnimngs
-        transparency = 0.0
-#endif
+
+        
     END SUBROUTINE usm_raytrace
     
@@ -3936 +3939 @@
         REAL(wp)                              :: acoef              !< actual coefficient of diurnal profile of anthropogenic heat
 
-#if defined( __mpi3 )
+
         dxdir = (/dz,dy,dy,dx,dx/)
         
@@ -4147 +4150 @@
        ENDIF
 
-#endif
+
     END SUBROUTINE usm_surface_energy_balance