Ignore:
Timestamp:
Apr 19, 2017 9:34:46 AM (4 years ago)
Author:
kanani
Message:

small bugfix, formatting and new PCM output

File:
1 edited

Legend:

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

    r2114 r2209  
    2121! Current revisions:
    2222! ------------------
    23 !
     23! cpp switch __mpi3 removed,
     24! minor formatting,
     25! small bugfix for division by zero (Krc)
    2426!
    2527! Former revisions:
     
    14891491        CALL location_message( '    calculation of SVF and CSF', .TRUE. )
    14901492
    1491 #if defined( __mpi3 )
     1493
    14921494!--     precalculate face areas for different face directions using normal vector
    14931495        DO d = 0, 9
     
    19641966        CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
    19651967       
    1966 #endif
     1968
    19671969    END SUBROUTINE usm_calc_svf
    19681970
     
    24872489!------------------------------------------------------------------------------!
    24882490    PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)
    2489         IMPLICIT NONE
    2490         REAL(wp), DIMENSION(3), INTENT(in)      :: origin    !< ray origin
    2491         REAL(wp), DIMENSION(3), INTENT(in)      :: uvect     !< ray unit vector
    2492         INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross  !< found boundary crossing (d, z, y, x)
    2493         REAL(wp), DIMENSION(3)                  :: crossdist !< crossing distance
    2494         INTEGER(iwp), DIMENSION(3)              :: bdyd      !< boundary direction
    2495         REAL(wp)                                :: bdydim    !<
    2496         REAL(wp)                                :: dist      !<
    2497         INTEGER(iwp)                            :: seldim    !< found fist crossing index
    2498         INTEGER(iwp)                            :: d         !<
    2499 
    2500         bdydim = nzut + .5_wp !< top boundary
    2501         bdyd(1) = isky
    2502         crossdist(1) = (bdydim - origin(1)) / uvect(1)
    2503 
    2504         IF ( uvect(2) >= 0._wp )  THEN
    2505             bdydim = ny + .5_wp !< north global boundary
    2506             bdyd(2) = inorthb
    2507         ELSE
    2508             bdydim = -.5_wp !< south global boundary
    2509             bdyd(2) = isouthb
    2510         ENDIF
    2511         crossdist(2) = (bdydim - origin(2)) / uvect(2)
    2512 
    2513         IF ( uvect(3) >= 0._wp )  THEN
    2514             bdydim = nx + .5_wp !< east global boundary
    2515             bdyd(3) = ieastb
    2516         ELSE
    2517             bdydim = -.5_wp !< west global boundary
    2518             bdyd(3) = iwestb
    2519         ENDIF
    2520         crossdist(3) = (bdydim - origin(3)) / uvect(3)
    2521 
    2522         seldim = minloc(crossdist, 1)
    2523         dist = crossdist(seldim)
    2524         d = bdyd(seldim)
    2525 
    2526         bdycross(1) = d
    2527         bdycross(2:4) = NINT( origin(:) + uvect(:)*dist &
    2528                         + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
     2491   
     2492       IMPLICIT NONE
     2493       
     2494       INTEGER(iwp) ::  d       !<
     2495       INTEGER(iwp) ::  seldim  !< found fist crossing index
     2496
     2497       INTEGER(iwp), DIMENSION(3)              ::  bdyd      !< boundary direction       
     2498       INTEGER(iwp), DIMENSION(4), INTENT(out) ::  bdycross  !< found boundary crossing (d, z, y, x)
     2499       
     2500       REAL(wp)                                ::  bdydim  !<
     2501       REAL(wp)                                ::  dist    !<
     2502       
     2503       REAL(wp), DIMENSION(3)             ::  crossdist  !< crossing distance
     2504       REAL(wp), DIMENSION(3), INTENT(in) ::  origin     !< ray origin
     2505       REAL(wp), DIMENSION(3), INTENT(in) ::  uvect      !< ray unit vector
     2506 
     2507
     2508       bdydim       = nzut + .5_wp  !< top boundary
     2509       bdyd(1)      = isky
     2510       crossdist(1) = ( bdydim - origin(1) ) / uvect(1)  !< subroutine called only when uvect(1)>0
     2511
     2512       IF ( uvect(2) == 0._wp )  THEN
     2513          crossdist(2) = huge(1._wp)
     2514       ELSE
     2515          IF ( uvect(2) >= 0._wp )  THEN
     2516             bdydim  = ny + .5_wp  !< north global boundary
     2517             bdyd(2) = inorthb
     2518          ELSE
     2519             bdydim  = -.5_wp  !< south global boundary
     2520             bdyd(2) = isouthb
     2521          ENDIF
     2522          crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
     2523       ENDIF
     2524
     2525       IF ( uvect(3) == 0._wp )  THEN
     2526          crossdist(3) = huge(1._wp)
     2527       ELSE
     2528          IF ( uvect(3) >= 0._wp )  THEN
     2529             bdydim  = nx + .5_wp  !< east global boundary
     2530             bdyd(3) = ieastb
     2531          ELSE
     2532             bdydim  = -.5_wp  !< west global boundary
     2533             bdyd(3) = iwestb
     2534          ENDIF
     2535          crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
     2536       ENDIF
     2537
     2538       seldim = minloc(crossdist, 1)
     2539       dist   = crossdist(seldim)
     2540       d      = bdyd(seldim)
     2541
     2542       bdycross(1)   = d
     2543       bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
     2544                                       + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
     2545                       
    25292546    END SUBROUTINE
    25302547
     
    28642881       urban_surface = .TRUE.
    28652882       
    2866 !
    2867 !--    Check whether pre-processor (cpp) option "__mpi3" is set. It is required
    2868 !--    for the full functionality of the USM. "__mpi3" directive is implemented,
    2869 !--    because some compilers cannot handle MPI-3 operations, hence, these parts
    2870 !--    of code shall only be compiled if explicitly enabled.
    2871 #if ! defined ( __mpi3 )
    2872           message_string = 'urban surface model requires compilation of ' //   &
    2873                            'PALM with pre-processor directive -D__mpi3'
    2874           CALL message( 'usm_parin', 'PA0503', 1, 2, 0, 6, 0 )
    2875 #endif
    2876 
    28772883
    28782884 10    CONTINUE
     
    32203226        REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
    32213227
    3222 #if defined( __mpi3 )
     3228
    32233229!--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
    32243230!--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
     
    33633369       
    33643370        visible = .TRUE.
    3365        
    3366 #else
    3367         visible      = .FALSE.                          !Set variables to avoid compiler warnimngs
    3368         transparency = 0.0
    3369 #endif
     3371
     3372       
    33703373    END SUBROUTINE usm_raytrace
    33713374   
     
    39363939        REAL(wp)                              :: acoef              !< actual coefficient of diurnal profile of anthropogenic heat
    39373940
    3938 #if defined( __mpi3 )
     3941
    39393942        dxdir = (/dz,dy,dy,dx,dx/)
    39403943       
     
    41474150       ENDIF
    41484151
    4149 #endif
     4152
    41504153    END SUBROUTINE usm_surface_energy_balance
    41514154
Note: See TracChangeset for help on using the changeset viewer.