Changeset 2024


Ignore:
Timestamp:
Oct 12, 2016 4:42:37 PM (5 years ago)
Author:
kanani
Message:

changes related to urban surface model and output of ssws

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r2012 r2024  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added missing CASE, error message and unit for ssws*,
     23! increased number of possible output quantities to 500.
    2324!
    2425! Former revisions:
     
    28342835    IF ( data_output_user(1) /= ' ' )  THEN
    28352836       i = 1
    2836        DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
     2837       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
    28372838          i = i + 1
    28382839       ENDDO
    28392840       j = 1
    2840        DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
    2841           IF ( i > 100 )  THEN
     2841       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
     2842          IF ( i > 500 )  THEN
    28422843             message_string = 'number of output quantitities given by data' // &
    2843                 '_output and data_output_user exceeds the limit of 100'
     2844                '_output and data_output_user exceeds the limit of 500'
    28442845             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
    28452846          ENDIF
     
    30183019             CONTINUE
    30193020
    3020           CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 't*', &
     3021          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 'ssws*', 't*', &
    30213022                 'u*', 'z0*', 'z0h*', 'z0q*' )
    30223023             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     
    30543055                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
    30553056             ENDIF
     3057             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
     3058                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3059                                 'res passive_scalar = .TRUE.'
     3060                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
     3061             ENDIF             
    30563062
    30573063             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
     
    30613067             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
    30623068             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
     3069             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'     
    30633070             IF ( TRIM( var ) == 't*'     )  unit = 'K'
    30643071             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r2012 r2024  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added missing lad_s initialization
    2323!
    2424! Former revisions:
     
    707707        REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
    708708
     709        lad_s = 0.0_wp
    709710        OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
    710711                action='READ', status='OLD', form='FORMATTED', err=515)
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r2012 r2024  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added missing CASE for ssws*
    2323!
    2424! Former revisions:
     
    343343                ENDIF
    344344                shf_av = 0.0_wp
     345               
     346             CASE ( 'ssws*' )
     347                IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
     348                   ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
     349                ENDIF
     350                ssws_av = 0.0_wp               
    345351
    346352             CASE ( 't*' )
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r2012 r2024  
    2121! Current revisions:
    2222! ------------------
    23 !
     23! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
     24! optimization of MPI-RMA operations,
     25! declaration of pcbl as integer,
     26! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
     27! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
     28! use of new control parameter varnamelength,
     29! added output variables usm_rad_ressw, usm_rad_reslw,
     30! minor formatting changes,
     31! minor optimizations.
    2432!
    2533! Former revisions:
     
    4351! ------------
    4452! 2016/6/9 - Initial version of the USM (Urban Surface Model)
    45 !            authors: Jaroslav Resler, Pavel Krc (CTU in Prague, ICS AS in Prague)
     53!            authors: Jaroslav Resler, Pavel Krc
     54!                     (Czech Technical University in Prague and Institute of
     55!                      Computer Science of the Czech Academy of Sciences, Prague)
    4656!            with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
    4757!            partly inspired by PALM LSM (B. Maronga)
     
    7383!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
    7484!>
     85!> @todo Check optimizations for RMA operations
     86!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
     87!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
     88!>       factor 3 between min and max time
    7589!------------------------------------------------------------------------------!
    7690 MODULE urban_surface_mod
     
    92106               time_since_reference_point, surface_pressure,                   &
    93107               g, pt_surface, large_scale_forcing, lsf_surf,                   &
    94                time_do3d, dt_do3d, average_count_3d, urban_surface
     108               time_do3d, dt_do3d, average_count_3d, varnamelength,            &
     109               urban_surface
    95110
    96111    USE cpulog,                                                                &
     
    161176    INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers in urban surface layer above top of buildings
    162177    INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
     178    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
    163179    INTEGER(iwp), PARAMETER                        ::  ndcsf = 2                          !< number of dimensions of real values in CSF
    164     INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF
     180    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
     181    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
    165182    INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
    166183    INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
     
    205222    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf
    206223                                                                        !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
    207     INTEGER(iwp)                                   ::  nsvfl            !< number of svf (excluding csf) for local processor
    208     INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor 
     224    INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
     225    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
    209226                                                                        !< needed only during calc_svf but must be here because it is
    210                                                                         !< shared between subroutines usm_calc_svf and usm_raytrace                                                                         
    211     INTEGER(iwp)                                   ::  nsvfcsfl         !< sum of svf+csf for local processor
    212 
     227                                                                        !< shared between subroutines usm_calc_svf and usm_raytrace
    213228
    214229!-- type for calculation of svf
     
    240255    INTEGER(iwp)                                   ::  msvf, mcsf       !< mod for swapping the growing array
    241256    INTEGER(iwp), PARAMETER                        ::  gasize = 10000   !< initial size of growing arrays
     257!-- temporary arrays for calculation of csf in raytracing
     258    INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
     259    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
     260    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
     261    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored
     262    INTEGER(kind=MPI_ADDRESS_KIND), &
     263                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
     264    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
    242265
    243266!-- arrays storing the values of USM
    244267    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]
    245     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors
    246                                                                         !< for individual local surfaces and plant canopy sinks
     268    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
    247269    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
    248270    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
     
    276298    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    277299    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     300    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
     301    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    278302    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface   
    279303   
    280304!-- block variables needed for calculation of the plant canopy model inside the urban surface model
    281     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  pcbl             !< z,y,x coordinates of i-th local plant canopy box pcbl[:,i] = [z, y, x]
    282     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[z,y,x]
     305    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
     306    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
     307                                                                        !< for local surfaces
     308    INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
     309    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[k,j,i]
    283310    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
    284311    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
     
    486513        INTEGER(iwp)                            :: i, j, k, d, l, ir, jr, ids
    487514        INTEGER(iwp)                            :: nzubl, nzutl, isurf, ipcgb
    488         INTEGER                                 :: procid
     515        INTEGER(iwp)                            :: procid
    489516
    490517       
     
    828855 
    829856        INTEGER(iwp)                                       :: i, j, k, l, ids, iwl,istat
    830         CHARACTER (len=20)                                 :: var, surfid
     857        CHARACTER (len=varnamelength)                      :: var, surfid
    831858        INTEGER(iwp), PARAMETER                            :: nd = 5
    832859        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
     
    861888           SELECT CASE ( TRIM( var ) )
    862889               
    863                 CASE ( 'usm_radnet' )
     890                CASE ( 'usm_rad_net' )
    864891!--                 array of complete radiation balance
    865892                    IF ( .NOT.  ALLOCATED(rad_net_av) )  THEN
     
    931958                    ENDIF
    932959
     960                CASE ( 'usm_rad_ressw' )
     961!--                 array of residua of sw radiation absorbed in surface after last reflection
     962                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
     963                        ALLOCATE( surfins_av(startenergy:endenergy) )
     964                        surfins_av = 0.0_wp
     965                    ENDIF
     966                                   
     967                CASE ( 'usm_rad_reslw' )
     968!--                 array of residua of lw radiation absorbed in surface after last reflection
     969                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
     970                        ALLOCATE( surfinl_av(startenergy:endenergy) )
     971                        surfinl_av = 0.0_wp
     972                    ENDIF
     973                                   
    933974                CASE ( 'usm_rad_hf' )
    934975!--                 array of heat flux from radiation for surfaces after i-th reflection
     
    9591000                        t_surf_av = 0.0_wp
    9601001                    ENDIF
    961                    
     1002
    9621003                CASE ( 'usm_t_wall' )
    9631004!--                 wall temperature for iwl layer of walls and land
     
    9761017           SELECT CASE ( TRIM( var ) )
    9771018               
    978                 CASE ( 'usm_radnet' )
     1019                CASE ( 'usm_rad_net' )
    9791020!--                 array of complete radiation balance
    9801021                    DO l = startenergy, endenergy
     
    10581099                    ENDDO
    10591100                   
     1101                CASE ( 'usm_rad_ressw' )
     1102!--                 array of residua of sw radiation absorbed in surface after last reflection
     1103                    DO l = startenergy, endenergy
     1104                        IF ( surfl(id,l) == ids )  THEN
     1105                            surfins_av(l) = surfins_av(l) + surfins(l)
     1106                        ENDIF
     1107                    ENDDO
     1108                                   
     1109                CASE ( 'usm_rad_reslw' )
     1110!--                 array of residua of lw radiation absorbed in surface after last reflection
     1111                    DO l = startenergy, endenergy
     1112                        IF ( surfl(id,l) == ids )  THEN
     1113                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
     1114                        ENDIF
     1115                    ENDDO
     1116                   
    10601117                CASE ( 'usm_rad_hf' )
    10611118!--                 array of heat flux from radiation for surfaces after i-th reflection
     
    11071164           SELECT CASE ( TRIM( var ) )
    11081165               
    1109                 CASE ( 'usm_radnet' )
     1166                CASE ( 'usm_rad_net' )
    11101167!--                 array of complete radiation balance
    11111168                    DO l = startenergy, endenergy
     
    11871244                    ENDDO
    11881245
     1246                CASE ( 'usm_rad_ressw' )
     1247!--                 array of residua of sw radiation absorbed in surface after last reflection
     1248                    DO l = startenergy, endenergy
     1249                        IF ( surfl(id,l) == ids )  THEN
     1250                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
     1251                        ENDIF
     1252                    ENDDO
     1253                                   
     1254                CASE ( 'usm_rad_reslw' )
     1255!--                 array of residua of lw radiation absorbed in surface after last reflection
     1256                    DO l = startenergy, endenergy
     1257                        IF ( surfl(id,l) == ids )  THEN
     1258                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
     1259                        ENDIF
     1260                    ENDDO
     1261                   
    11891262                CASE ( 'usm_rad_hf' )
    11901263!--                 array of heat flux from radiation for surfaces after i-th reflection
     
    12181291                        ENDIF
    12191292                    ENDDO
    1220                    
     1293
    12211294                CASE ( 'usm_t_wall' )
    12221295!--                 wall temperature for  iwl layer of walls and land
     
    12261299                        ENDIF
    12271300                    ENDDO
    1228                
     1301
    12291302           END SELECT
    12301303
     
    13911464        REAL(wp), DIMENSION(3)                      :: uv
    13921465        LOGICAL                                     :: visible
    1393         REAL(wp)                                    :: sz, sy, sx, tz, ty, tx, transparency, rirrf, sqdist, svfsum
    1394         !REAL(wp)                                    :: rsvf
     1466        REAL(wp), DIMENSION(3)                      :: sa, ta          !< real coordinates z,y,x of source and target
     1467        REAL(wp)                                    :: transparency, rirrf, sqdist, svfsum
    13951468        INTEGER(iwp)                                :: isurflt, isurfs, isurflt_prev
    13961469        INTEGER(iwp)                                :: itx, ity, itz
     
    14391512        IF ( plant_canopy )  THEN
    14401513            ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
     1514            maxboxesg = nx + ny + nzu + 1
     1515!--         temporary arrays storing values for csf calculation during raytracing
     1516            ALLOCATE( boxes(3, maxboxesg) )
     1517            ALLOCATE( crlens(maxboxesg) )
     1518
    14411519#if defined( __parallel )
    14421520            ALLOCATE( planthl(nys:nyn,nxl:nxr) )
     
    14451523            CALL MPI_AllGather( planthl, nnx*nny, MPI_INTEGER, &
    14461524                                plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
    1447             DEALLOCATE(planthl)
     1525            DEALLOCATE( planthl )
     1526           
     1527!--         temporary arrays storing values for csf calculation during raytracing
     1528            ALLOCATE( lad_ip(maxboxesg) )
     1529            ALLOCATE( lad_disp(maxboxesg) )
    14481530
    14491531            IF ( usm_lad_rma )  THEN
     1532                ALLOCATE( lad_s_ray(maxboxesg) )
     1533               
     1534                ! set conditions for RMA communication
    14501535                CALL MPI_Info_create(minfo, ierr)
    14511536                CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
     
    15001585                IF ( i < nxl  .OR.  i > nxr &
    15011586                     .OR.  j < nys  .OR.  j > nyn ) CYCLE
    1502                 tx = REAL(i)
    1503                 ty = REAL(j)
    1504                 tz = REAL(k)
     1587                ta = (/ REAL(k), REAL(j), REAL(i) /)
    15051588
    15061589                DO isurfs = 1, nsurf
     
    15121595                     
    15131596                    sd = surf(id, isurfs)
    1514                     sz = REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd)
    1515                     sy = REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd)
    1516                     sx = REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)
     1597                    sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
     1598                            REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
     1599                            REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
    15171600
    15181601!--                 unit vector source -> target
    1519                     uv = (/ (tz-sz)*dz, (ty-sy)*dy, (tx-sx)*dx /)
     1602                    uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
    15201603                    sqdist = SUM(uv(:)**2)
    15211604                    uv = uv / SQRT(sqdist)
     
    15271610
    15281611!--                 raytrace while not creating any canopy sink factors
    1529                     CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, 1._wp, .FALSE., &
     1612                    CALL usm_raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., &
    15301613                            visible, transparency, win_lad)
    15311614                    IF ( .NOT.  visible ) CYCLE
     
    15581641            td = surfl(id, isurflt)
    15591642            IF ( td >= isky  .AND.  .NOT.  plant_canopy ) CYCLE
    1560             tz = REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td)
    1561             ty = REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td)
    1562             tx = REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td)
     1643            ta = (/ REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td),  &
     1644                      REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td),  &
     1645                      REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td)  /)
    15631646            DO isurfs = 1, nsurf
    15641647                IF ( .NOT.  usm_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
     
    15701653                 
    15711654                sd = surf(id, isurfs)
    1572                 sz = REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd)
    1573                 sy = REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd)
    1574                 sx = REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)
     1655                sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
     1656                        REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
     1657                        REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
    15751658
    15761659!--             unit vector source -> target
    1577                 uv = (/ (tz-sz)*dz, (ty-sy)*dy, (tx-sx)*dx /)
     1660                uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
    15781661                sqdist = SUM(uv(:)**2)
    15791662                uv = uv / SQRT(sqdist)
     
    15861669
    15871670!--             raytrace + process plant canopy sinks within
    1588                 CALL usm_raytrace((/sz,sy,sx/), (/tz,ty,tx/), isurfs, rirrf, facearea(td), .TRUE., &
     1671                CALL usm_raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
    15891672                        visible, transparency, win_lad)
    15901673               
     
    16231706
    16241707        CALL location_message( '    waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
    1625 
     1708!--     deallocate temporary global arrays
     1709        DEALLOCATE(nzterr)
     1710       
     1711        IF ( plant_canopy )  THEN
     1712!--         finalize mpi_rma communication and deallocate temporary arrays
    16261713#if defined( __parallel )
    1627         IF ( plant_canopy )  THEN
    16281714            IF ( usm_lad_rma )  THEN
    16291715                CALL MPI_Win_flush_all(win_lad, ierr)
     
    16321718!--             free MPI window
    16331719                CALL MPI_Win_free(win_lad, ierr)
     1720               
     1721!--             deallocate temporary arrays storing values for csf calculation during raytracing
     1722                DEALLOCATE( lad_s_ray )
     1723!--             usm_lad is the pointer to lad_s_rma in case of usm_lad_rma
     1724!--             and must not be deallocated here
    16341725            ELSE
    16351726                DEALLOCATE(usm_lad)
    16361727                DEALLOCATE(usm_lad_g)
    16371728            ENDIF
     1729#else
     1730            DEALLOCATE(usm_lad)
     1731#endif
     1732            DEALLOCATE( boxes )
     1733            DEALLOCATE( crlens )
     1734            DEALLOCATE( plantt )
    16381735        ENDIF
    1639 #else
    1640         DEALLOCATE(usm_lad)
    1641 #endif
    1642 
    1643 !--     deallocate temporary global arrays
    1644         IF ( ALLOCATED(nzterr) ) DEALLOCATE(nzterr)
    1645         IF ( ALLOCATED(plantt) ) DEALLOCATE(plantt)
    1646        
     1736
     1737        CALL location_message( '    calculation of the complete SVF array', .TRUE. )
     1738
    16471739!--     sort svf ( a version of quicksort )
    16481740        CALL quicksort_svf(asvf,1,nsvfl)
    1649        
     1741
     1742        ALLOCATE( svf(ndsvf,nsvfl) )
     1743        ALLOCATE( svfsurf(idsvf,nsvfl) )
     1744
     1745        !< load svf from the structure array to plain arrays
     1746        isurflt_prev = -1
     1747        ksvf = 1
     1748        svfsum = 0._wp
     1749        DO isvf = 1, nsvfl
     1750!--         normalize svf per target face
     1751            IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
     1752                IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
     1753!--                 TODO detect and log when normalization differs too much from 1
     1754                    svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
     1755                ENDIF
     1756                isurflt_prev = asvf(ksvf)%isurflt
     1757                isvf_surflt = isvf
     1758                svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
     1759            ELSE
     1760                svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
     1761            ENDIF
     1762
     1763            svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
     1764            svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
     1765
     1766!--         next element
     1767            ksvf = ksvf + 1
     1768        ENDDO
     1769
     1770        IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
     1771!--         TODO detect and log when normalization differs too much from 1
     1772            svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
     1773        ENDIF
     1774
     1775!--     deallocate temporary asvf array
     1776!--     DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target
     1777!--     via pointing pointer - we need to test original targets
     1778        IF ( ALLOCATED(asvf1) )  THEN
     1779            DEALLOCATE(asvf1)
     1780        ENDIF
     1781        IF ( ALLOCATED(asvf2) )  THEN
     1782            DEALLOCATE(asvf2)
     1783        ENDIF
     1784
    16501785        npcsfl = 0
    16511786        IF ( plant_canopy )  THEN
     1787
     1788            CALL location_message( '    calculation of the complete CSF array', .TRUE. )
     1789
    16521790!--         sort and merge csf for the last time, keeping the array size to minimum
    16531791            CALL usm_merge_and_grow_csf(-1)
     
    17201858!--         scatter and gather the number of elements to and from all processor
    17211859!--         and calculate displacements
    1722             CALL mpi_alltoall(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
     1860            CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
    17231861           
    17241862            npcsfl = SUM(ipcsflt)
     
    17881926                npcsfl = kcsf
    17891927            ENDIF
    1790            
    1791         ENDIF  !< plant_canopy
    1792        
    1793         CALL location_message( '    calculation of the complete SVF array', .TRUE. )
    1794         nsvfcsfl = nsvfl + npcsfl
    1795        
    1796         ALLOCATE( svf(ndsvf,nsvfcsfl) )
    1797         ALLOCATE( svfsurf(ndsvf,nsvfcsfl) )
    1798 
    1799         !< load svf from the structure array to plain arrays
    1800         isurflt_prev = -1
    1801         ksvf = 1
    1802         svfsum = 0._wp
    1803         DO isvf = 1, nsvfl
    1804 !--         normalize svf per target face
    1805             IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
    1806                 IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    1807 !--                 TODO detect and log when normalization differs too much from 1
    1808                     svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
    1809                 ENDIF
    1810                 isurflt_prev = asvf(ksvf)%isurflt
    1811                 isvf_surflt = isvf
    1812                 svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
    1813             ELSE
    1814                 svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
    1815             ENDIF
    1816            
    1817             svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
    1818             svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
    1819                    
    1820 !--         next element
    1821             ksvf = ksvf + 1   
    1822         ENDDO
    1823        
    1824         IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    1825 !--         TODO detect and log when normalization differs too much from 1
    1826             svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
    1827         ENDIF
    1828        
    1829 !--     deallocate temporary asvf array
    1830 !--     DEALLOCATE(asvf) - ifort has a problem with deallocation of allocatable target
    1831 !--     via pointing pointer - we need to test original targets
    1832         IF ( ALLOCATED(asvf1) )  THEN
    1833             DEALLOCATE(asvf1)
    1834         ENDIF
    1835         IF ( ALLOCATED(asvf2) )  THEN
    1836             DEALLOCATE(asvf2)
    1837         ENDIF
    1838 
    1839         IF ( plant_canopy )  THEN
    1840             CALL location_message( '    calculation of the complete CSF part of the array', .TRUE. )
    1841             IF ( npcsfl > 0 )  THEN
    1842                 DO isvf = 1, npcsfl
    1843                     svf(:,nsvfl+isvf) = pcsflt(:,isvf)
    1844                     svfsurf(1,nsvfl+isvf) =  gridpcbl(kpcsflt(1,isvf),kpcsflt(2,isvf),kpcsflt(3,isvf))
    1845                     svfsurf(2,nsvfl+isvf) =  kpcsflt(4,isvf)
     1928
     1929            ncsfl = npcsfl
     1930            IF ( ncsfl > 0 )  THEN
     1931                ALLOCATE( csf(ndcsf,ncsfl) )
     1932                ALLOCATE( csfsurf(idcsf,ncsfl) )
     1933                DO icsf = 1, ncsfl
     1934                    csf(:,icsf) = pcsflt(:,icsf)
     1935                    csfsurf(1,icsf) =  gridpcbl(kpcsflt(1,icsf),kpcsflt(2,icsf),kpcsflt(3,icsf))
     1936                    csfsurf(2,icsf) =  kpcsflt(4,icsf)
    18461937                ENDDO
    18471938            ENDIF
     
    18771968        CHARACTER (len=*),INTENT(OUT)   ::  unit     !:
    18781969       
    1879         CHARACTER (len=20)              :: var
     1970        CHARACTER (len=varnamelength)   :: var
    18801971
    18811972        var = TRIM(variable)
    1882         IF ( var(1:11)=='usm_radnet_'  .OR.  var(1:13)=='usm_rad_insw_'  .OR.           &
    1883              var(1:13)=='usm_rad_inlw_'  .OR.  var(1:16)=='usm_rad_inswdir_'  .OR.      &
    1884              var(1:16)=='usm_rad_inswdif_'  .OR.  var(1:16)=='usm_rad_inswref_'  .OR.   &
    1885              var(1:16)=='usm_rad_inlwdif_'  .OR.  var(1:16)=='usm_rad_inlwref_'  .OR.   &
    1886              var(1:14)=='usm_rad_outsw_'  .OR.  var(1:14)=='usm_rad_outlw_'  .OR.       &
    1887              var(1:11)=='usm_rad_hf_'  .OR.                                             &
    1888              var(1:9) =='usm_wshf_'  .OR.  var(1:9)=='usm_wghf_' )  THEN
     1973        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.        &
     1974             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.    &
     1975             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
     1976             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
     1977             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.     &
     1978             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.     &
     1979             var(1:11) == 'usm_rad_hf_'  .OR.                                             &
     1980             var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' )  THEN
    18891981            unit = 'W/m2'
    1890         ELSE IF ( var(1:10) =='usm_t_surf'  .OR.  var(1:10) =='usm_t_wall' )  THEN
     1982        ELSE IF ( var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall' )  THEN
    18911983            unit = 'K'
    1892         ELSE IF ( var(1:9) =='usm_surfz'  .OR.  var(1:7) =='usm_svf'  .OR.              &
    1893                   var(1:7) =='usm_dif'  .OR.  var(1:11) =='usm_surfcat'  .OR.           &
    1894                   var(1:11) =='usm_surfalb'  .OR.  var(1:12) =='usm_surfemis')  THEN
     1984        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              &
     1985                  var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
     1986                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis')  THEN
    18951987            unit = '1'
    1896         ELSE IF ( plant_canopy  .AND.  var(1:7) =='usm_lad' )  THEN
     1988        ELSE IF ( plant_canopy  .AND.  var(1:7) == 'usm_lad' )  THEN
    18971989            unit = 'm2/m3'
    1898         ELSE IF ( plant_canopy  .AND.  var(1:14) == 'usm_canopy_khf' )  THEN
     1990        ELSE IF ( plant_canopy  .AND.  var(1:13) == 'usm_canopy_hr' )  THEN
    18991991            unit = 'K/s'
    19001992        ELSE
     
    19722064        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)     ::  temp_pf    !< temp array for urban surface output procedure
    19732065       
    1974         CHARACTER (len=20)                                     :: var, surfid
     2066        CHARACTER (len=varnamelength)                          :: var, surfid
    19752067        INTEGER(iwp), PARAMETER                                :: nd = 5
    19762068        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
     
    20852177              ENDDO
    20862178
    2087           CASE ( 'usm_radnet' )
     2179          CASE ( 'usm_rad_net' )
    20882180!--           array of complete radiation balance
    20892181              DO isurf = dirstart(ids), dirend(ids)
     
    22022294                   ELSE
    22032295                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
     2296                   ENDIF
     2297                 ENDIF
     2298              ENDDO
     2299
     2300          CASE ( 'usm_rad_ressw' )
     2301!--           average of array of residua of sw radiation absorbed in surface after last reflection
     2302              DO isurf = dirstart(ids), dirend(ids)
     2303                 IF ( surfl(id,isurf) == ids )  THEN
     2304                   IF ( av == 0 )  THEN
     2305                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
     2306                   ELSE
     2307                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
     2308                   ENDIF
     2309                 ENDIF
     2310              ENDDO
     2311
     2312          CASE ( 'usm_rad_reslw' )
     2313!--           average of array of residua of lw radiation absorbed in surface after last reflection
     2314              DO isurf = dirstart(ids), dirend(ids)
     2315                 IF ( surfl(id,isurf) == ids )  THEN
     2316                   IF ( av == 0 )  THEN
     2317                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
     2318                   ELSE
     2319                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
    22042320                   ENDIF
    22052321                 ENDIF
     
    22772393              ENDDO
    22782394             
    2279           CASE ( 'usm_canopy_khf' )
    2280 !--           canopy kinematic heat flux
     2395          CASE ( 'usm_canopy_hr' )
     2396!--           canopy heating rate
    22812397              DO i = nxl, nxr
    22822398                 DO j = nys, nyn
     
    23222438        CHARACTER (len=*), INTENT(OUT) ::  grid_z      !<
    23232439
    2324         CHARACTER (len=20)            :: var
     2440        CHARACTER (len=varnamelength)  :: var
    23252441
    23262442        var = TRIM(variable)
    2327         IF ( var(1:11)=='usm_radnet_'  .OR.  var(1:13) =='usm_rad_insw_'  .OR.            &
    2328              var(1:13) =='usm_rad_inlw_'  .OR.  var(1:16) =='usm_rad_inswdir_'  .OR.      &
    2329              var(1:16) =='usm_rad_inswdif_'  .OR.  var(1:16) =='usm_rad_inswref_'  .OR.   &
    2330              var(1:16) =='usm_rad_inlwdif_'  .OR.  var(1:16) =='usm_rad_inlwref_'  .OR.   &
    2331              var(1:14) =='usm_rad_outsw_'  .OR.  var(1:14) =='usm_rad_outlw_'  .OR.       &
    2332              var(1:11) =='usm_rad_hf_'  .OR.                                              &
    2333              var(1:9) == 'usm_wshf_'  .OR.  var(1:9)== 'usm_wghf_'  .OR.                  &
    2334              var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall'  .OR.             &
    2335              var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                   &
    2336              var(1:7) =='usm_dif'  .OR.  var(1:11) =='usm_surfcat'  .OR.                  &
    2337              var(1:11) =='usm_surfalb'  .OR.  var(1:12) =='usm_surfemis'  .OR.            &
    2338              var(1:7) == 'usm_lad'  .OR.  var(1:14) == 'usm_canopy_khf' )  THEN
     2443        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.          &
     2444             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.      &
     2445             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR.   &
     2446             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR.   &
     2447             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
     2448             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
     2449             var(1:11) == 'usm_rad_hf_'  .OR.                                               &
     2450             var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
     2451             var(1:10) == 'usm_t_surf'  .OR.  var(1:10) == 'usm_t_wall'  .OR.               &
     2452             var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     &
     2453             var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
     2454             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
     2455             var(1:7) == 'usm_lad'  .OR.  var(1:13) == 'usm_canopy_hr' )  THEN
    23392456
    23402457            found = .TRUE.
     
    25312648       
    25322649        IF ( read_svf_on_init )  THEN
    2533 !--         read svf and svfsurf data from file
     2650!--         read svf, csf, svfsurf and csfsurf data from file
    25342651            CALL location_message( '    Start reading SVF from file', .TRUE. )
    25352652            CALL usm_read_svf_from_file()
     
    25452662
    25462663        IF ( write_svf_on_init )  THEN
    2547 !--         write svf and svfsurf data to file
     2664!--         write svf, csf svfsurf and csfsurf data to file
     2665            CALL location_message( '    Store SVF and CSF to file', .TRUE. )
    25482666            CALL usm_write_svf_to_file()
    25492667        ENDIF
     
    26862804! Description:
    26872805! ------------
     2806!> Parin for &usm_par for urban surface model
     2807!------------------------------------------------------------------------------!
     2808    SUBROUTINE usm_parin
     2809
     2810       IMPLICIT NONE
     2811
     2812       CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
     2813
     2814       NAMELIST /urban_surface_par/                                            &
     2815                           land_category,                                      &
     2816                           mrt_factors,                                        &
     2817                           nrefsteps,                                          &
     2818                           pedestrant_category,                                &
     2819                           ra_horiz_coef,                                      &
     2820                           read_svf_on_init,                                   &
     2821                           roof_category,                                      &
     2822                           split_diffusion_radiation,                          &
     2823                           urban_surface,                                      &
     2824                           usm_anthropogenic_heat,                             &
     2825                           usm_energy_balance_land,                            &
     2826                           usm_energy_balance_wall,                            &
     2827                           usm_material_model,                                 &
     2828                           usm_lad_rma,                                        &
     2829                           wall_category,                                      &
     2830                           write_svf_on_init
     2831
     2832       line = ' '
     2833
     2834!
     2835!--    Try to find urban surface model package
     2836       REWIND ( 11 )
     2837       line = ' '
     2838       DO   WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
     2839          READ ( 11, '(A)', END=10 )  line
     2840       ENDDO
     2841       BACKSPACE ( 11 )
     2842
     2843!
     2844!--    Read user-defined namelist
     2845       READ ( 11, urban_surface_par )
     2846
     2847!
     2848!--    Set flag that indicates that the land surface model is switched on
     2849       urban_surface = .TRUE.
     2850
     2851
     2852 10    CONTINUE
     2853
     2854    END SUBROUTINE usm_parin
     2855
     2856
     2857!------------------------------------------------------------------------------!
     2858! Description:
     2859! ------------
    26882860!> This subroutine calculates interaction of the solar radiation
    26892861!> with urban surface and updates surface, roofs and walls heatfluxes.
     
    26952867       
    26962868        INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep
    2697         INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, ipcgb
     2869        INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
    26982870        INTEGER(iwp), DIMENSION(4) :: bdycross
    26992871        REAL(wp), DIMENSION(3,3)   :: mrot            !< grid rotation matrix (xyz)
     
    28363008!--         pcsf first pass
    28373009            isurf1 = -1  !< previous processed pcgb
    2838             DO isvf = nsvfl+1, nsvfcsfl
    2839                 ipcgb = svfsurf(1, isvf)
     3010            DO icsf = 1, ncsfl
     3011                ipcgb = csfsurf(1, icsf)
    28403012                i = pcbl(ix,ipcgb)
    28413013                j = pcbl(iy,ipcgb)
    28423014                k = pcbl(iz,ipcgb)
    2843                 isurfsrc = svfsurf(2, isvf)
     3015                isurfsrc = csfsurf(2, icsf)
    28443016
    28453017                IF ( zenith(0) > 0  .AND.  ipcgb /= isurf1 )  THEN
     
    29183090
    29193091!--         radiation absorbed by plant canopy
    2920             DO isvf = nsvfl+1, nsvfcsfl
    2921                 ipcgb = svfsurf(1, isvf)
    2922                 isurfsrc = svfsurf(2, isvf)
     3092            DO icsf = 1, ncsfl
     3093                ipcgb = csfsurf(1, icsf)
     3094                isurfsrc = csfsurf(2, icsf)
    29233095
    29243096                IF ( surf(id, isurfsrc) < isky )  THEN
    2925                     pcbinsw(ipcgb) = pcbinsw(ipcgb) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
    2926 !--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc)
     3097                    pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
     3098!--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc)
    29273099                ENDIF
    29283100            ENDDO
     
    29983170        REAL(wp), INTENT(out)                  :: transparency !< along whole path
    29993171        INTEGER(iwp), INTENT(in)               :: win_lad
    3000         INTEGER(iwp)                           :: k, d
     3172        INTEGER(iwp)                           :: i, j, k, d
    30013173        INTEGER(iwp)                           :: seldim       !< dimension to be incremented
    30023174        INTEGER(iwp)                           :: ncsb         !< no of written plant canopy sinkboxes
     
    30173189        INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before
    30183190                                                               !< the processor in the question
    3019         INTEGER(iwp)                           :: ig, ip
    3020         REAL(wp)                               :: lad_s_target
    3021         INTEGER(kind=MPI_ADDRESS_KIND)         :: lad_disp
    3022         REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp
     3191        INTEGER(iwp)                           :: ip           !< number of processor where gridbox reside
     3192        INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
     3193        REAL(wp)                               :: lad_s_target !< recieved lad_s of particular grid box
     3194        REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
    30233195
    30243196!--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
     
    30863258                IF ( box(1) <= nzterr(ig) )  THEN
    30873259                    visible = .FALSE.
    3088                     IF ( ncsb > 0 )  THEN
    3089 !--                     rewind written plant canopy sink factors - they are invalid
    3090                         ncsfl = ncsfl - ncsb
    3091                     ENDIF
    30923260                    RETURN
    30933261                ENDIF
     
    30953263                IF ( plant_canopy )  THEN
    30963264                    IF ( box(1) <= plantt(ig) )  THEN
     3265                        ncsb = ncsb + 1
     3266                        boxes(:,ncsb) = box
     3267                        crlens(ncsb) = crlen
    30973268#if defined( __parallel )
    3098                         lad_disp = (box(3)-px*nnx)*(nny*nzu) + (box(2)-py*nny)*nzu + box(1)-nzub
    3099                         IF ( usm_lad_rma )  THEN
    3100 !--                         Read LAD using MPI RMA
    3101                             CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
    3102                             CALL MPI_Get(lad_s_target, 1, MPI_REAL, ip, lad_disp, 1, MPI_REAL, &
    3103                                 win_lad, ierr)
    3104                             IF ( ierr /= 0 )  THEN
    3105                                 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
    3106                                 CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
    3107                             ENDIF
    3108                            
    3109                             CALL MPI_Win_flush_local(ip, win_lad, ierr)
    3110                             IF ( ierr /= 0 )  THEN
    3111                                 WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local'
    3112                                 CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
    3113                             ENDIF
    3114                             CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
    3115                         ELSE
    3116                             lad_s_target = usm_lad_g(ip*nnx*nny*nzu + lad_disp)
    3117                         ENDIF
    3118 #else
    3119                         lad_s_target = usm_lad(box(1),box(2),box(3))
     3269                        lad_ip(ncsb) = ip
     3270                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzu) + (box(2)-py*nny)*nzu + box(1)-nzub
    31203271#endif
    3121                         cursink = 1._wp - exp(-ext_coef * lad_s_target &
    3122                             * crlen*realdist)
    3123 
    3124                         IF ( create_csf )  THEN
    3125 !--                         write svf values into the array
    3126                             ncsb = ncsb + 1
    3127                             ncsfl = ncsfl + 1
    3128                             acsf(ncsfl)%ip = ip
    3129                             acsf(ncsfl)%itx = box(3)
    3130                             acsf(ncsfl)%ity = box(2)
    3131                             acsf(ncsfl)%itz = box(1)
    3132                             acsf(ncsfl)%isurfs = isrc
    3133                             acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
    3134                             acsf(ncsfl)%rtransp = REAL(transparency, wp)
    3135                         ENDIF  !< create_csf
    3136 
    3137                         transparency = transparency * (1._wp - cursink)
    3138                        
    31393272                    ENDIF
    31403273                ENDIF
     
    31473280        ENDDO
    31483281       
     3282        IF ( plant_canopy )  THEN
     3283#if defined( __parallel )
     3284            IF ( usm_lad_rma )  THEN
     3285!--             send requests for lad_s to appropriate processor
     3286                CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
     3287                DO i = 1, ncsb
     3288                    CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
     3289                                 1, MPI_REAL, win_lad, ierr)
     3290                    IF ( ierr /= 0 )  THEN
     3291                        WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
     3292                        CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
     3293                    ENDIF
     3294                ENDDO
     3295               
     3296!--             wait for all pending local requests complete
     3297                CALL MPI_Win_flush_local_all(win_lad, ierr)
     3298                IF ( ierr /= 0 )  THEN
     3299                    WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
     3300                    CALL message( 'usm_raytrace', 'PA0519', 1, 2, 0, 6, 0 )
     3301                ENDIF
     3302                CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
     3303               
     3304            ENDIF
     3305#endif
     3306
     3307!--         calculate csf and transparency
     3308            DO i = 1, ncsb
     3309#if defined( __parallel )
     3310                IF ( usm_lad_rma )  THEN
     3311                    lad_s_target = lad_s_ray(i)
     3312                ELSE
     3313                    lad_s_target = usm_lad_g(lad_ip(i)*nnx*nny*nzu + lad_disp(i))
     3314                ENDIF
     3315#else
     3316                lad_s_target = usm_lad(boxes(1,i),boxes(2,i),boxes(3,i))
     3317#endif
     3318                cursink = 1._wp - exp(-ext_coef * lad_s_target * crlens(i)*realdist)
     3319
     3320                IF ( create_csf )  THEN
     3321!--                 write svf values into the array
     3322                    ncsfl = ncsfl + 1
     3323                    acsf(ncsfl)%ip = lad_ip(i)
     3324                    acsf(ncsfl)%itx = boxes(3,i)
     3325                    acsf(ncsfl)%ity = boxes(2,i)
     3326                    acsf(ncsfl)%itz = boxes(1,i)
     3327                    acsf(ncsfl)%isurfs = isrc
     3328                    acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
     3329                    acsf(ncsfl)%rtransp = REAL(transparency, wp)
     3330                ENDIF  !< create_csf
     3331
     3332                transparency = transparency * (1._wp - cursink)
     3333               
     3334            ENDDO
     3335        ENDIF
     3336       
    31493337        visible = .TRUE.
     3338       
    31503339    END SUBROUTINE usm_raytrace
    31513340   
     
    31973386           
    31983387#if defined( __parallel ) && ! defined ( __check )
    3199             CALL mpi_barrier( comm2d, ierr )
     3388            CALL MPI_BARRIER( comm2d, ierr )
    32003389#endif
    32013390        ENDDO
     
    32313420           
    32323421#if defined( __parallel ) && ! defined ( __check )
    3233             CALL mpi_barrier( comm2d, ierr )
     3422            CALL MPI_BARRIER( comm2d, ierr )
    32343423#endif
    32353424        ENDDO
     
    32563445       CHARACTER (LEN=30) ::  variable_chr  !< dummy variable to read string
    32573446       
    3258        INTEGER            ::  i             !< running index     
    3259        
     3447       INTEGER(iwp)       ::  i             !< running index
     3448
    32603449
    32613450       DO  i = 0, io_blocks-1
     
    33173506
    33183507        IMPLICIT NONE
    3319         INTEGER                      :: fsvf = 89
    3320         INTEGER                      :: i
     3508        INTEGER(iwp)                 :: fsvf = 89
     3509        INTEGER(iwp)                 :: i
    33213510        CHARACTER(usm_version_len)   :: usm_version_field
    33223511        CHARACTER(svf_code_len)      :: svf_code_field
     
    33363525                ENDIF
    33373526               
    3338 !--             read nsvfcsfl, nsvfl
    3339                 READ ( fsvf ) nsvfcsfl, nsvfl
    3340                 IF ( nsvfcsfl <= 0 )  THEN
     3527!--             read nsvfl, ncsfl
     3528                READ ( fsvf ) nsvfl, ncsfl
     3529                IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
    33413530                    WRITE( message_string, * ) 'Wrong number of SVF or CSF'
    33423531                    CALL message( 'usm_read_svf_from_file', 'UI0012', 1, 2, 0, 6, 0 )
    33433532                ELSE
    3344                     WRITE(message_string,*) '    Number of SVF and CSF to read', nsvfcsfl, nsvfl
     3533                    WRITE(message_string,*) '    Number of SVF and CSF to read', nsvfl, ncsfl
    33453534                    CALL location_message( message_string, .TRUE. )
    33463535                ENDIF
    33473536               
    3348                 ALLOCATE(svf(ndsvf,nsvfcsfl))
    3349                 ALLOCATE(svfsurf(ndsvf,nsvfcsfl))
    3350                
     3537                ALLOCATE(svf(ndsvf,nsvfl))
     3538                ALLOCATE(svfsurf(idsvf,nsvfl))
    33513539                READ(fsvf) svf
    33523540                READ(fsvf) svfsurf
     3541                IF ( plant_canopy )  THEN
     3542                    ALLOCATE(csf(ndcsf,ncsfl))
     3543                    ALLOCATE(csfsurf(idcsf,ncsfl))
     3544                    READ(fsvf) csf
     3545                    READ(fsvf) csfsurf
     3546                ENDIF
    33533547                READ ( fsvf ) svf_code_field
    33543548               
     
    33623556            ENDIF
    33633557#if defined( __parallel )
    3364             CALL mpi_barrier( comm2d, ierr )
     3558            CALL MPI_BARRIER( comm2d, ierr )
    33653559#endif
    33663560        ENDDO
     
    33793573   
    33803574        CHARACTER(12)                                         :: wtn
    3381         INTEGER                                               :: wtc
     3575        INTEGER(iwp)                                          :: wtc
    33823576        REAL(wp), DIMENSION(n_surface_params)                 :: wtp
    33833577   
     
    34903684            ENDIF
    34913685#if defined( __parallel ) && ! defined ( __check )
    3492             CALL mpi_barrier( comm2d, ierr )
     3686            CALL MPI_BARRIER( comm2d, ierr )
    34933687#endif
    34943688        ENDDO
     
    39124106       THEN
    39134107#if defined( __parallel )
    3914           IF ( collective_wait )  CALL mpi_barrier( comm2d, ierr )
     4108          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    39154109              CALL mpi_allreduce( force_radiation_call_l, force_radiation_call,         &
    39164110                                  1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
     
    39344128       IMPLICIT NONE
    39354129
    3936        INTEGER, INTENT(IN) :: mod_count
    3937        INTEGER :: i
     4130       INTEGER(iwp), INTENT(IN) :: mod_count
     4131       INTEGER(iwp)            :: i
    39384132     
    39394133#if defined( __nopointer )
     
    40274221       IMPLICIT NONE
    40284222       
    4029        INTEGER ::  i
     4223       INTEGER(iwp) ::  i
    40304224
    40314225       DO  i = 0, io_blocks-1
     
    40584252! Description:
    40594253! ------------
    4060 !> Subroutine stores svf and svfsurf data to files
     4254!> Subroutine stores svf, svfsurf, csf and csfsurf data to a file.
    40614255!------------------------------------------------------------------------------!
    40624256    SUBROUTINE usm_write_svf_to_file
    4063    
     4257
    40644258        IMPLICIT NONE
    4065         INTEGER             :: fsvf = 89
    4066         INTEGER             :: i
    4067        
     4259        INTEGER(iwp)        :: fsvf = 89
     4260        INTEGER(iwp)        :: i
     4261
    40684262        DO  i = 0, io_blocks-1
    40694263            IF ( i == io_group )  THEN
     
    40724266
    40734267                WRITE ( fsvf )  usm_version
    4074                 WRITE ( fsvf )  nsvfcsfl, nsvfl
     4268                WRITE ( fsvf )  nsvfl, ncsfl
    40754269                WRITE ( fsvf )  svf
    40764270                WRITE ( fsvf )  svfsurf
     4271                IF ( plant_canopy )  THEN
     4272                    WRITE ( fsvf )  csf
     4273                    WRITE ( fsvf )  csfsurf
     4274                ENDIF
    40774275                WRITE ( fsvf )  TRIM(svf_code)
    4078                
     4276
    40794277                CLOSE (fsvf)
    40804278#if defined( __parallel )
    4081                 CALL mpi_barrier( comm2d, ierr )
     4279                CALL MPI_BARRIER( comm2d, ierr )
    40824280#endif
    40834281            ENDIF
     
    40864284
    40874285
    4088 !------------------------------------------------------------------------------!
    4089 ! Description:
    4090 ! ------------
    4091 !> Parin for &usm_par for urban surface model
    4092 !------------------------------------------------------------------------------!
    4093     SUBROUTINE usm_parin
    4094 
    4095        IMPLICIT NONE
    4096 
    4097        CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
    4098 
    4099        NAMELIST /urban_surface_par/                                            & 
    4100                            land_category,                                      &
    4101                            mrt_factors,                                        &
    4102                            nrefsteps,                                          &
    4103                            pedestrant_category,                                &
    4104                            ra_horiz_coef,                                      &
    4105                            read_svf_on_init,                                   &
    4106                            roof_category,                                      &
    4107                            split_diffusion_radiation,                          &
    4108                            urban_surface,                                      &
    4109                            usm_anthropogenic_heat,                             &
    4110                            usm_energy_balance_land,                            &
    4111                            usm_energy_balance_wall,                            &
    4112                            usm_material_model,                                 &
    4113                            usm_lad_rma,                                        &
    4114                            wall_category,                                      &
    4115                            write_svf_on_init
    4116 
    4117        line = ' '
    4118        
    4119 !
    4120 !--    Try to find urban surface model package
    4121        REWIND ( 11 )
    4122        line = ' '
    4123        DO   WHILE ( INDEX( line, '&urban_surface_par' ) == 0 )
    4124           READ ( 11, '(A)', END=10 )  line
    4125        ENDDO
    4126        BACKSPACE ( 11 )
    4127 
    4128 !
    4129 !--    Read user-defined namelist
    4130        READ ( 11, urban_surface_par )
    4131 
    4132 !
    4133 !--    Set flag that indicates that the land surface model is switched on
    4134        urban_surface = .TRUE.
    4135        
    4136 
    4137  10    CONTINUE
    4138    
    4139     END SUBROUTINE usm_parin   
    4140 
    4141    
    41424286!------------------------------------------------------------------------------!
    41434287!
     
    41464290!> Block of auxiliary subroutines:
    41474291!> 1. quicksort and corresponding comparison
    4148 !> 2. usm_merge_and_grow_csf for implementation of "dynamical growing" array
    4149 !>    for svf and csf
     4292!> 2. usm_merge_and_grow_csf for implementation of "dynamical growing"
     4293!>    array for csf
    41504294!------------------------------------------------------------------------------!   
    41514295    PURE FUNCTION svf_lt(svf1,svf2) result (res)
     
    41684312        IMPLICIT NONE
    41694313        TYPE(t_svf), DIMENSION(:), INTENT(INOUT)  :: svfl
    4170         INTEGER, INTENT(IN)                       :: first, last
     4314        INTEGER(iwp), INTENT(IN)                  :: first, last
    41714315        TYPE(t_svf)                               :: x, t
    4172         INTEGER                                   :: i, j
     4316        INTEGER(iwp)                              :: i, j
    41734317
    41744318        IF ( first>=last ) RETURN
     
    42174361        IMPLICIT NONE
    42184362        TYPE(t_csf), DIMENSION(:), INTENT(INOUT)  :: csfl
    4219         INTEGER, INTENT(IN)                       :: first, last
     4363        INTEGER(iwp), INTENT(IN)                  :: first, last
    42204364        TYPE(t_csf)                               :: x, t
    4221         INTEGER                                   :: i, j
     4365        INTEGER(iwp)                              :: i, j
    42224366
    42234367        IF ( first>=last ) RETURN
     
    42434387   
    42444388    SUBROUTINE usm_merge_and_grow_csf(newsize)
    4245         INTEGER(iwp), INTENT(in)            :: newsize  !< new array size after grow, must be >= ncsfl
    4246                                                         !< or -1 to shrink to minimum
    4247         INTEGER(iwp)                        :: iread, iwrite
    4248         TYPE(t_csf), DIMENSION(:), POINTER  :: acsfnew
     4389        INTEGER(iwp), INTENT(in)                :: newsize  !< new array size after grow, must be >= ncsfl
     4390                                                            !< or -1 to shrink to minimum
     4391        INTEGER(iwp)                            :: iread, iwrite
     4392        TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
    42494393
    42504394        IF ( newsize == -1 )  THEN
     
    43254469        INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT)  :: kpcsflt
    43264470        REAL(wp), DIMENSION(:,:), INTENT(INOUT)      :: pcsflt
    4327         INTEGER, INTENT(IN)                          :: first, last
     4471        INTEGER(iwp), INTENT(IN)                     :: first, last
    43284472        REAL(wp), DIMENSION(ndcsf)                   :: t2
    43294473        INTEGER(iwp), DIMENSION(kdcsf)               :: x, t1
    4330         INTEGER                                      :: i, j
     4474        INTEGER(iwp)                                 :: i, j
    43314475
    43324476        IF ( first>=last ) RETURN
Note: See TracChangeset for help on using the changeset viewer.