Changeset 2977


Ignore:
Timestamp:
Apr 17, 2018 10:27:57 AM (6 years ago)
Author:
kanani
Message:

Fixes for radiative transfer model

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

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

    r2938 r2977  
    2525! -----------------
    2626! $Id$
     27! Implement changes from branch radiation (r2948-2971) with minor modifications
     28! (moh.hefny):
     29! - set radiation_interactions according to the existence of urban/land vertical
     30!   surfaces and trees to activiate RTM
     31! - set average_radiation to TRUE if RTM is activiated
     32!
     33! 2938 2018-03-27 15:52:42Z suehring
    2734! - Revise Inifor initialization for geostrophic wind components
    2835! - Initialize synthetic turbulence generator in case of Inifor initialization 
     
    506513
    507514    USE radiation_model_mod,                                                   &
    508         ONLY:  radiation_init, radiation, radiation_control, radiation_scheme, &
     515        ONLY:  average_radiation,                                              &
     516               radiation_init, radiation, radiation_control, radiation_scheme, &
    509517               radiation_calc_svf, radiation_write_svf,                        &
    510518               radiation_interaction, radiation_interactions,                  &
    511519               radiation_interaction_init, radiation_read_svf,                 &
    512                radiation_presimulate_solar_pos
     520               radiation_presimulate_solar_pos, radiation_interactions_on
    513521   
    514522    USE random_function_mod
     
    534542    USE surface_mod,                                                           &
    535543        ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
    536                 surf_usm_h, get_topography_top_index_ji
     544                surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist
    537545   
    538546    USE transpose_indices
     
    23572365!-- If required, set chemical emissions
    23582366!-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model)
    2359     IF ( air_chemistry ) THEN
     2367    IF ( air_chemistry )  THEN
    23602368       CALL chem_emissions
    23612369    ENDIF
     
    23642372!-- Initialize radiation processes
    23652373    IF ( radiation )  THEN
     2374!
     2375!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
     2376!--    The namelist parameter radiation_interactions_on can override this behavior.
     2377!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
     2378!--    init_surface_arrays.)
     2379       IF ( radiation_interactions_on )  THEN
     2380          IF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
     2381             radiation_interactions    = .TRUE.
     2382             average_radiation         = .TRUE.
     2383          ELSE
     2384             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
     2385                                                   !< calculations necessary in case of flat surface
     2386          ENDIF
     2387       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
     2388          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
     2389                           'vertical surfaces and/or trees exist. The model will run ' // &
     2390                           'without RTM (no shadows, no radiation reflections)'
     2391          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
     2392       ENDIF
    23662393!
    23672394!--    If required, initialize radiation interactions between surfaces
  • palm/trunk/SOURCE/palm.f90

    r2951 r2977  
    2525! -----------------
    2626! $Id$
     27! Deduct spinup_time from RUN_CONTROL output of main 3d run
     28! (use time_since_reference_point instead of simulated_time)
     29!
     30! 2951 2018-04-06 09:05:08Z kanani
    2731! Add log_point_s for pmci_init
    2832!
     
    240244               nest_domain, neutral, nudging, passive_scalar, runnr,           &
    241245               simulated_time, simulated_time_chr, spinup,                     &
     246               time_since_reference_point,                                     &
    242247               user_interface_current_revision,                                &
    243248               user_interface_required_revision, version, wall_heatflux,       &
     
    464469!
    465470!-- Set start time in format hh:mm:ss
    466     simulated_time_chr = time_to_string( simulated_time )
     471    simulated_time_chr = time_to_string( time_since_reference_point )
    467472
    468473!
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r2932 r2977  
    1919!
    2020! Current revisions:
    21 ! -----------------
     21! ------------------
    2222!
    2323!
     
    2525! -----------------
    2626! $Id$
     27! Implement changes from branch radiation (r2948-2971) with minor modifications,
     28! plus some formatting.
     29! (moh.hefny):
     30! Add plant canopy type to account for changes in LAD (based on the changes
     31! done by Resler & Pavel) and correct the error message to PALM Standard.
     32!
     33! 2932 2018-03-26 09:39:22Z maronga
    2734! renamed canopy_par to plant_canopy_parameters
    2835!
     
    213220    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
    214221    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
     222
     223    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
     224                                                !< of plant canopy (e.g. deciduous tree during winter)
    215225
    216226    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
     
    962972                                  alpha_lad, beta_lad, canopy_drag_coeff,      &
    963973                                  canopy_mode, cthf,                           &
    964                                   lad_surface,                                 &
     974                                  lad_surface, lad_type_coef,                  &
    965975                                  lad_vertical_gradient,                       &
    966976                                  lad_vertical_gradient_level,                 &
     
    971981       NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
    972982                                  canopy_mode, cthf,                           &
    973                                   lad_surface,                                 &
     983                                  lad_surface, lad_type_coef,                  &
    974984                                  lad_vertical_gradient,                       &
    975985                                  lad_vertical_gradient_level,                 &
     
    10341044!>
    10351045!> num_levels
    1036 !> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
    1037 !> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
    1038 !> dtype,x,y,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
     1046!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
     1047!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
     1048!> dtype,x,y,pctype,value(nzb),value(nzb+1), ... ,value(nzb+num_levels-1)
    10391049!> ...
    10401050!>
     
    10601070
    10611071       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
     1072       INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
    10621073       INTEGER(iwp)                        ::  i, j      !< running index
    10631074       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
     
    10731084!
    10741085!--    Open and read plant canopy input data
    1075        OPEN(152, file='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
    1076                  access='SEQUENTIAL', action='READ', status='OLD',             &
    1077                  form='FORMATTED', err=515)
    1078        READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
     1086       OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
     1087                 ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
     1088                 FORM='FORMATTED', ERR=515)
     1089       READ(152, *, ERR=516, END=517) nzp   !< read first line = number of vertical layers
    10791090       
    1080        ALLOCATE(col(0:nzp-1))
     1091       ALLOCATE( col(0:nzp-1) )
    10811092
    10821093       DO
    1083           READ(152, *, err=516, end=517) dtype, i, j, col(:)
    1084           IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
    1085 
    1086              SELECT CASE (dtype)
    1087                 CASE( 1 )   !< leaf area density
    1088 !
    1089 !--                 This is just the pure canopy layer assumed to be grounded to
    1090 !--                 a flat domain surface. At locations where plant canopy sits
    1091 !--                 on top of any kind of topography, the vertical plant column
    1092 !--                 must be "lifted", which is done in SUBROUTINE pcm_tendency.
    1093                     lad_s(0:nzp-1, j, i) = col(0:nzp-1)
    1094                    
    1095                 CASE DEFAULT
    1096                    WRITE(message_string, '(a,i2,a)')   &
    1097                       'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
    1098                    CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
    1099              END SELECT
     1094          READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
     1095          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
     1096         
     1097          SELECT CASE (dtype)
     1098             CASE( 1 )   !< leaf area density
     1099!
     1100!--             This is just the pure canopy layer assumed to be grounded to
     1101!--             a flat domain surface. At locations where plant canopy sits
     1102!--             on top of any kind of topography, the vertical plant column
     1103!--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
     1104                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
     1105                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
     1106                                              'Allowed values 0 <= pctype <= 10, ' //  &
     1107                                              'but pctype is ', pctype
     1108                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
     1109                ENDIF
     1110                lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype)
     1111               
     1112             CASE DEFAULT
     1113                WRITE(message_string, '(a,i2,a)')   &
     1114                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
     1115                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
     1116          END SELECT
    11001117       ENDDO
    11011118
     
    11071124
    11081125517    CLOSE(152)
    1109        DEALLOCATE(col)
     1126       DEALLOCATE( col )
    11101127       
    11111128       CALL exchange_horiz( lad_s, nbgp )
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2967 r2977  
    2222!
    2323! Current revisions:
    24 ! -----------------
     24! ------------------
    2525!
    2626!
     
    2828! -----------------
    2929! $Id$
     30! Implement changes from branch radiation (r2948-2971) with minor modifications,
     31! plus some formatting.
     32! (moh.hefny):
     33! - replaced plant_canopy by npcbl to check tree existence to avoid weird
     34!   allocation of related arrays (after domain decomposition some domains
     35!   contains no trees although plant_canopy (global parameter) is still TRUE).
     36! - added a namelist parameter to force RTM settings
     37! - enabled the option to switch radiation reflections off
     38! - renamed surf_reflections to surface_reflections
     39! - removed average_radiation flag from the namelist (now it is implicitly set
     40!   in init_3d_model according to RTM)
     41! - edited read and write sky view factors and CSF routines to account for
     42!   the sub-domains which may not contain any of them
     43!
     44! 2967 2018-04-13 11:22:08Z raasch
    3045! bugfix: missing parallel cpp-directives added
    3146!
     
    432447                sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
    433448                average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
    434                 radiation_interactions = .FALSE.,     & !< flag to control if radiation interactions via sky-view factors shall be considered
    435                 surf_reflections = .TRUE.              !< flag to switch the calculation of radiation interaction between surfaces.
     449                radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
     450                surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
    436451                                                        !< When it switched off, only the effect of buildings and trees shadow will
    437452                                                        !< will be considered. However fewer SVFs are expected.
     453                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
    438454
    439455
     
    680696    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
    681697    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
    682     INTEGER(iwp)                                   ::  npcbl            !< number of the plant canopy gridboxes in local processor
     698    INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
    683699    INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
    684700    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
     
    955971           surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,     &
    956972           surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
    957            pcbl, npcbl, iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
     973           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
    958974           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
    959975           nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
     
    961977           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    962978           radiation_interactions, startwall, startland, endland, endwall,     &
    963            skyvf, skyvft
     979           skyvf, skyvft, radiation_interactions_on, average_radiation
    964980
    965981#if defined ( __rrtmg )
     
    13091325             CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
    13101326          ENDIF
    1311        ENDIF
    1312 
    1313 !
    1314 !--    Radiation interactions
    1315        IF ( nrefsteps < 1  .AND.  radiation_interactions )  THEN
    1316           message_string = 'nrefsteps must be > 0 when using LSM/USM to' //    &
    1317                            'account for surface outgoing SW flux.'       //    &
    1318                            'You may set surf_reflections = .FALSE. to '  //    &
    1319                            'diable surface reflections instead.'
    1320           CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
    13211327       ENDIF
    13221328
     
    26342640                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
    26352641                                  dist_max_svf,                                &
    2636                                   average_radiation,                           &
    2637                                   surf_reflections, svfnorm_report_thresh
     2642                                  surface_reflections, svfnorm_report_thresh,  &
     2643                                  radiation_interactions_on
    26382644   
    26392645       NAMELIST /radiation_parameters/   albedo, albedo_type, albedo_lw_dir,   &
     
    26472653                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
    26482654                                  dist_max_svf,                                &
    2649                                   average_radiation,                           &
    2650                                   surf_reflections, svfnorm_report_thresh
     2655                                  surface_reflections, svfnorm_report_thresh,  &
     2656                                  radiation_interactions_on
    26512657   
    26522658       line = ' '
     
    26932699       radiation = .TRUE.
    26942700
     2701       IF ( .NOT.  radiation_interactions_on  .AND.  surface_reflections )  THEN
     2702          message_string = 'surface_reflections is allowed only when '      // &
     2703               'radiation_interactions_on is set to TRUE'
     2704          CALL message( 'radiation_parin', 'PA0487',1, 2, 0, 6, 0 )
     2705       ENDIF
     2706
    26952707 12    CONTINUE
    2696        
    2697 
    2698 
    2699 !--    Set radiation_interactions flag according to urban_ and land_surface flag
    2700        IF ( urban_surface  .OR.  land_surface )  radiation_interactions = .TRUE.
    27012708       
    27022709    END SUBROUTINE radiation_parin
     
    43544361         sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
    43554362
    4356          IF ( plant_canopy )  THEN
     4363         IF ( npcbl > 0 )  THEN
    43574364!--            precompute effective box depth with prototype Leaf Area Density
    43584365            pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
     
    44654472#endif
    44664473
    4467      DO isvf = 1, nsvfl
    4468          isurf = svfsurf(1, isvf)
    4469          k = surfl(iz, isurf)
    4470          j = surfl(iy, isurf)
    4471          i = surfl(ix, isurf)
    4472          isurfsrc = svfsurf(2, isvf)
    4473 
    4474 !--         for surface-to-surface factors we calculate thermal radiation in 1st pass
    4475          surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
    4476      ENDDO
     4474     IF ( surface_reflections)  THEN
     4475        DO  isvf = 1, nsvfl
     4476           isurf = svfsurf(1, isvf)
     4477           k     = surfl(iz, isurf)
     4478           j     = surfl(iy, isurf)
     4479           i     = surfl(ix, isurf)
     4480           isurfsrc = svfsurf(2, isvf)
     4481!
     4482!--        For surface-to-surface factors we calculate thermal radiation in 1st pass
     4483           surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
     4484        ENDDO
     4485     ENDIF
    44774486
    44784487     !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors
     
    44944503     ENDIF
    44954504
    4496      IF ( plant_canopy )  THEN
     4505     IF ( npcbl > 0 )  THEN
    44974506
    44984507         pcbinswdir(:) = 0._wp
     
    45364545!        surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
    45374546
     4547     IF ( .NOT.  surface_reflections )  THEN
     4548!
     4549!--     Set nrefsteps to 0 to disable reflections       
     4550        nrefsteps = 0
     4551        surfoutsl = albedo_surf * surfins
     4552        surfoutll = (1._wp - emiss_surf) * surfinl
     4553        surfoutsw = surfoutsw + surfoutsl
     4554        surfoutlw = surfoutlw + surfoutll
     4555     ENDIF
     4556
    45384557!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    45394558!--     Next passes - reflections
     
    45854604
    45864605!--  push heat flux absorbed by plant canopy to respective 3D arrays
    4587      IF ( plant_canopy )  THEN
     4606     IF ( npcbl > 0 )  THEN
    45884607         pc_heating_rate(:,:,:) = 0._wp
    45894608         DO ipcgb = 1, npcbl
     
    47334752!--  domain when using average_radiation in the respective radiation model
    47344753
    4735      IF ( average_radiation )  THEN
    4736 
    4737 !--     precalculate face areas for all face directions using normal vector
    4738         DO d = 0, nsurf_type
    4739             facearea(d) = 1._wp
    4740             IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
    4741             IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
    4742             IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
    4743         ENDDO
    4744 !
    4745 !--     absorbed/received SW & LW and emitted LW energy of all physical
    4746 !--     surfaces (land and urban) in local processor
    4747         pinswl = 0._wp
    4748         pinlwl = 0._wp
    4749         pabsswl = 0._wp
    4750         pabslwl = 0._wp
    4751         pemitlwl = 0._wp
    4752         emiss_sum_surfl = 0._wp
    4753         area_surfl = 0._wp
    4754         DO  i = 1, nsurfl
    4755            d = surfl(id, i)
    4756 !--        received SW & LW
    4757            pinswl = pinswl + surfinsw(i) * facearea(d)
    4758            pinlwl = pinlwl + surfinlw(i) * facearea(d)
    4759 !--        absorbed SW & LW
    4760            pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
    4761                                                   surfinsw(i) * facearea(d)
    4762            pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
    4763 !--        emitted LW
    4764            pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
    4765 !--        emissivity and area sum
    4766            emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
    4767            area_surfl = area_surfl + facearea(d)
    4768         END DO
    4769 !
    4770 !--     add the absorbed SW energy by plant canopy
    4771         IF ( plant_canopy )  THEN
    4772            pabsswl = pabsswl + SUM(pcbinsw)
    4773            pabslwl = pabslwl + SUM(pcbinlw)
    4774         ENDIF
    4775 !
    4776 !--     gather all rad flux energy in all processors
     4754!--  Precalculate face areas for all face directions using normal vector
     4755     DO d = 0, nsurf_type
     4756        facearea(d) = 1._wp
     4757        IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
     4758        IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
     4759        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
     4760     ENDDO
     4761!
     4762!--  absorbed/received SW & LW and emitted LW energy of all physical
     4763!--  surfaces (land and urban) in local processor
     4764     pinswl = 0._wp
     4765     pinlwl = 0._wp
     4766     pabsswl = 0._wp
     4767     pabslwl = 0._wp
     4768     pemitlwl = 0._wp
     4769     emiss_sum_surfl = 0._wp
     4770     area_surfl = 0._wp
     4771     DO  i = 1, nsurfl
     4772        d = surfl(id, i)
     4773!--  received SW & LW
     4774        pinswl = pinswl + surfinsw(i) * facearea(d)
     4775        pinlwl = pinlwl + surfinlw(i) * facearea(d)
     4776!--   absorbed SW & LW
     4777        pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
     4778                                                surfinsw(i) * facearea(d)
     4779        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
     4780!--   emitted LW
     4781        pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
     4782!--   emissivity and area sum
     4783        emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
     4784        area_surfl = area_surfl + facearea(d)
     4785     END DO
     4786!
     4787!--  add the absorbed SW energy by plant canopy
     4788     IF ( npcbl > 0 )  THEN
     4789        pabsswl = pabsswl + SUM(pcbinsw)
     4790        pabslwl = pabslwl + SUM(pcbinlw)
     4791     ENDIF
     4792!
     4793!--  gather all rad flux energy in all processors
    47774794#if defined( __parallel )
    4778         CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
    4779         CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
    4780         CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    4781         CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    4782         CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    4783         CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    4784         CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     4795     CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     4796     CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     4797     CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     4798     CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     4799     CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     4800     CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     4801     CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    47854802#else
    4786         pinsw = pinswl
    4787         pinlw = pinlwl
    4788         pabssw = pabsswl
    4789         pabslwl = pabslw
    4790         pemitlwl = pemitlw
    4791         emiss_sum_surf = emiss_sum_surfl
    4792         area_surf = area_surfl
     4803     pinsw = pinswl
     4804     pinlw = pinlwl
     4805     pabssw = pabsswl
     4806     pabslwl = pabslw
     4807     pemitlwl = pemitlw
     4808     emiss_sum_surf = emiss_sum_surfl
     4809     area_surf = area_surfl
    47934810#endif
    47944811
    4795 !--     (1) albedo
    4796         IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
    4797 
    4798 !--     (2) average emmsivity
    4799         IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
    4800 
    4801 !--     (3) temperature
    4802         t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
    4803 
    4804      ENDIF
     4812!--  (1) albedo
     4813     IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
     4814
     4815!--  (2) average emmsivity
     4816     IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
     4817
     4818!--  (3) temperature
     4819     t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
     4820     
    48054821
    48064822    CONTAINS
     
    50715087
    50725088!--    fill gridpcbl and pcbl
    5073        IF ( plant_canopy )  THEN
     5089       IF ( npcbl > 0 )  THEN
    50745090           ALLOCATE( pcbl(iz:ix, 1:npcbl) )
    50755091           ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) )
     
    55935609
    55945610            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp )
    5595 
    5596             DO isurfs = 1, nsurf
    5597                 IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
    5598                     surfl(iz, isurflt), surfl(id, isurflt), &
    5599                     surf(ix, isurfs), surf(iy, isurfs), &
    5600                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
    5601                     CYCLE
    5602                 ENDIF
     5611!
     5612!--         Following calculations only required for surface_reflections
     5613            IF ( surface_reflections )  THEN
     5614
     5615               DO  isurfs = 1, nsurf
     5616                  IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
     5617                     surfl(iz, isurflt), surfl(id, isurflt), &
     5618                     surf(ix, isurfs), surf(iy, isurfs), &
     5619                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
     5620                     CYCLE
     5621                  ENDIF
    56035622                 
    5604                 sd = surf(id, isurfs)
    5605                 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
    5606                         REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
    5607                         REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
    5608 
    5609 !--             unit vector source -> target
    5610                 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
    5611                 sqdist = SUM(uv(:)**2)
    5612                 uv = uv / SQRT(sqdist)
    5613 
    5614 !--             reject raytracing above max distance
    5615                 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
    5616                     ray_skip_maxdist = ray_skip_maxdist + 1
    5617                     CYCLE
    5618                 ENDIF
    5619                
    5620 !--             irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
    5621                 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
    5622                     * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
    5623                     / (pi * sqdist) & ! square of distance between centers
    5624                     * facearea(sd)
    5625 
    5626 !--             reject raytracing for potentially too small view factor values
    5627                 IF ( rirrf < min_irrf_value ) THEN
    5628                     ray_skip_minval = ray_skip_minval + 1
    5629                     CYCLE
    5630                 ENDIF
    5631 
    5632 !--             raytrace + process plant canopy sinks within
    5633                 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
    5634                         visible, transparency, win_lad)
    5635 
    5636                 IF ( .NOT.  visible ) CYCLE
    5637                 ! rsvf = rirrf * transparency
    5638 
    5639 !--             write to the svf array
    5640                 nsvfl = nsvfl + 1
    5641 !--             check dimmension of asvf array and enlarge it if needed
    5642                 IF ( nsvfla < nsvfl )  THEN
    5643                     k = nsvfla * 2
    5644                     IF ( msvf == 0 )  THEN
     5623                  sd = surf(id, isurfs)
     5624                  sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
     5625                          REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
     5626                          REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
     5627
     5628!--               unit vector source -> target
     5629                  uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
     5630                  sqdist = SUM(uv(:)**2)
     5631                  uv = uv / SQRT(sqdist)
     5632
     5633!--               reject raytracing above max distance
     5634                  IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
     5635                     ray_skip_maxdist = ray_skip_maxdist + 1
     5636                     CYCLE
     5637                  ENDIF
     5638                 
     5639!--               irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
     5640                  rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
     5641                      * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
     5642                      / (pi * sqdist) & ! square of distance between centers
     5643                      * facearea(sd)
     5644
     5645!--               reject raytracing for potentially too small view factor values
     5646                  IF ( rirrf < min_irrf_value ) THEN
     5647                      ray_skip_minval = ray_skip_minval + 1
     5648                      CYCLE
     5649                  ENDIF
     5650
     5651!--               raytrace + process plant canopy sinks within
     5652                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
     5653                                visible, transparency, win_lad)
     5654
     5655                  IF ( .NOT.  visible ) CYCLE
     5656                 ! rsvf = rirrf * transparency
     5657
     5658!--               write to the svf array
     5659                  nsvfl = nsvfl + 1
     5660!--               check dimmension of asvf array and enlarge it if needed
     5661                  IF ( nsvfla < nsvfl )  THEN
     5662                     k = nsvfla * 2
     5663                     IF ( msvf == 0 )  THEN
    56455664                        msvf = 1
    56465665                        ALLOCATE( asvf1(k) )
     
    56485667                        asvf1(1:nsvfla) = asvf2
    56495668                        DEALLOCATE( asvf2 )
    5650                     ELSE
     5669                     ELSE
    56515670                        msvf = 0
    56525671                        ALLOCATE( asvf2(k) )
     
    56545673                        asvf2(1:nsvfla) = asvf1
    56555674                        DEALLOCATE( asvf1 )
    5656                     ENDIF
    5657 
    5658 !                     WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
    5659 !                     CALL radiation_write_debug_log( msg )
    5660                    
    5661                     nsvfla = k
    5662                 ENDIF
    5663 !--             write svf values into the array
    5664                 asvf(nsvfl)%isurflt = isurflt
    5665                 asvf(nsvfl)%isurfs = isurfs
    5666                 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency
    5667                 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor
    5668             ENDDO
     5675                     ENDIF
     5676
     5677!                      WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
     5678!                      CALL radiation_write_debug_log( msg )
     5679                     
     5680                     nsvfla = k
     5681                  ENDIF
     5682!--               write svf values into the array
     5683                  asvf(nsvfl)%isurflt = isurflt
     5684                  asvf(nsvfl)%isurfs = isurfs
     5685                  asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency
     5686                  asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor
     5687               ENDDO
     5688            ENDIF
    56695689        ENDDO
    56705690
     
    67336753 !--             read nsvfl, ncsfl
    67346754              READ ( fsvf ) nsvfl, ncsfl, nsurfl_from_file
    6735               IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
     6755              IF ( nsvfl < 0  .OR.  ncsfl < 0 )  THEN
    67366756                  WRITE( message_string, * ) 'Wrong number of SVF or CSF'
    67376757                  CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
     
    67476767                  CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
    67486768              ENDIF
    6749              
    6750               IF ( .NOT. ALLOCATED( skyvf ) )    ALLOCATE( skyvf(nsurfl) )
    6751               IF ( .NOT. ALLOCATED( skyvft ) )   ALLOCATE( skyvft(nsurfl) )
    6752               IF ( .NOT. ALLOCATED( svf ) )      ALLOCATE( svf(ndsvf,nsvfl) )
    6753               IF ( .NOT. ALLOCATED( svfsurf ) )  ALLOCATE( svfsurf(idsvf,nsvfl) )
    6754               READ(fsvf) skyvf
    6755               READ(fsvf) skyvft
    6756               READ(fsvf) svf
    6757               READ(fsvf) svfsurf
    6758               IF ( plant_canopy )  THEN
    6759                   IF ( .NOT. ALLOCATED( csf ) )      ALLOCATE( csf(ndcsf,ncsfl) )
    6760                   IF ( .NOT. ALLOCATED( csfsurf ) )  ALLOCATE( csfsurf(idcsf,ncsfl) )
    6761                   READ(fsvf) csf
    6762                   READ(fsvf) csfsurf
     6769              IF ( nsurfl > 0 )  THEN
     6770                 IF ( .NOT. ALLOCATED( skyvf ) )    ALLOCATE( skyvf(nsurfl) )
     6771                 IF ( .NOT. ALLOCATED( skyvft ) )   ALLOCATE( skyvft(nsurfl) )
     6772                 READ(fsvf) skyvf
     6773                 READ(fsvf) skyvft
     6774              ENDIF
     6775               
     6776              IF ( nsvfl > 0 )  THEN
     6777                 IF ( .NOT. ALLOCATED( svf ) )      ALLOCATE( svf(ndsvf,nsvfl) )
     6778                 IF ( .NOT. ALLOCATED( svfsurf ) )  ALLOCATE( svfsurf(idsvf,nsvfl) )
     6779                 READ(fsvf) svf
     6780                 READ(fsvf) svfsurf
     6781              ENDIF
     6782
     6783              IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
     6784                 IF ( .NOT. ALLOCATED( csf ) )      ALLOCATE( csf(ndcsf,ncsfl) )
     6785                 IF ( .NOT. ALLOCATED( csfsurf ) )  ALLOCATE( csfsurf(idcsf,ncsfl) )
     6786                 READ(fsvf) csf
     6787                 READ(fsvf) csfsurf
    67636788              ENDIF
    67646789              READ ( fsvf ) svf_code_field
    6765              
     6790               
    67666791              IF ( TRIM(svf_code_field) /= TRIM(svf_code) )  THEN
    6767                   WRITE( message_string, * ) 'Wrong structure of binary svf file'
    6768                   CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
     6792                 WRITE( message_string, * ) 'Wrong structure of binary svf file'
     6793                 CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
    67696794              ENDIF       
    67706795!
    6771 !--          Close binary file                 
    6772              CALL close_file( fsvf )
     6796!--           Close binary file                 
     6797              CALL close_file( fsvf )
    67736798               
    67746799           ENDIF
    67756800#if defined( __parallel )
    6776             CALL MPI_BARRIER( comm2d, ierr )
     6801           CALL MPI_BARRIER( comm2d, ierr )
    67776802#endif
    67786803        ENDDO
     
    68016826                WRITE ( fsvf )  rad_version
    68026827                WRITE ( fsvf )  nsvfl, ncsfl, nsurfl
    6803                 WRITE ( fsvf )  skyvf
    6804                 WRITE ( fsvf )  skyvft
    6805                 WRITE ( fsvf )  svf
    6806                 WRITE ( fsvf )  svfsurf
    6807                 IF ( plant_canopy )  THEN
     6828                IF ( nsurfl > 0 ) THEN
     6829                   WRITE ( fsvf )  skyvf
     6830                   WRITE ( fsvf )  skyvft
     6831                ENDIF
     6832                IF ( nsvfl > 0 ) THEN
     6833                   WRITE ( fsvf )  svf
     6834                   WRITE ( fsvf )  svfsurf
     6835                ENDIF
     6836                IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
    68086837                    WRITE ( fsvf )  csf
    68096838                    WRITE ( fsvf )  csfsurf
  • palm/trunk/SOURCE/surface_mod.f90

    r2970 r2977  
    2626! -----------------
    2727! $Id$
     28! Implement changes from branch radiation (r2948-2971) with minor modifications,
     29! plus some formatting.
     30! (moh.hefny)
     31! Added flag to check the existence of vertical urban/land surfaces, required
     32! to activate RTM
     33!
     34! 2970 2018-04-13 15:09:23Z suehring
    2835! Remove un-necessary initialization of surface elements in old large-scale
    2936! forcing mode
     
    450457    INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data
    451458
     459    LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces
     460                                                    !< in the domain (required to activiate RTM)
     461
    452462
    453463    SAVE
     
    488498    PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file,       &
    489499           ns_v_on_file, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,       &
    490            surf_usm_h, surf_usm_v, surf_type
     500           surf_usm_h, surf_usm_v, surf_type, vertical_surfaces_exist
    491501!
    492502!-- Public subroutines and functions
     
    615625    SUBROUTINE init_surface_arrays
    616626
     627
     628       USE pegrid
     629
     630
    617631       IMPLICIT NONE
    618632
     
    628642       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
    629643       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces
     644
     645       INTEGER(iwp)              ::  num_surf_v_l !< number of vertically-aligned local urban/land surfaces
     646       INTEGER(iwp)              ::  num_surf_v   !< number of vertically-aligned total urban/land surfaces
    630647
    631648       LOGICAL ::  building                       !< flag indicating building grid point
     
    920937                                               nys, nyn, nxl, nxr )
    921938       ENDDO
    922 
     939!
     940!--    Set the flag for the existence of vertical urban/land surfaces
     941       num_surf_v_l = 0
     942       DO  l = 0, 3
     943          num_surf_v_l = num_surf_v_l + surf_usm_v(l)%ns + surf_lsm_v(l)%ns
     944       ENDDO
     945
     946#if defined( __parallel )
     947       CALL MPI_ALLREDUCE( num_surf_v_l, num_surf_v, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr)
     948#else
     949       num_surf_v = num_surf_v_l
     950#endif
     951       IF ( num_surf_v > 0 ) vertical_surfaces_exist = .TRUE.
     952       
    923953
    924954    END SUBROUTINE init_surface_arrays
  • palm/trunk/SOURCE/time_integration.f90

    r2941 r2977  
    2525! -----------------
    2626! $Id$
     27! Implement changes from branch radiation (r2948-2971) with minor modifications.
     28! (moh.hefny)
     29! Fixed bug in if statement
     30!
     31! 2941 2018-04-03 11:54:58Z kanani
    2732! Deduct spinup_time from RUN_CONTROL output of main 3d run
    2833! (use time_since_reference_point instead of simulated_time)
     
    10341039                CALL cpu_log( log_point(50), 'radiation', 'stop' )
    10351040
    1036                 IF ( urban_surface  .OR.  land_surface  .AND.                  &
     1041                IF ( ( urban_surface  .OR.  land_surface )  .AND.               &
    10371042                     radiation_interactions )  THEN
    10381043                   CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r2963 r2977  
    2828! -----------------
    2929! $Id$
     30! Implement changes from branch radiation (r2948-2971) with minor modifications.
     31! (moh.hefny):
     32! Extended exn for all model domain height to avoid the need to get nzut.
     33!
     34! 2963 2018-04-12 14:47:44Z suehring
    3035! Introduce index for vegetation/wall, pavement/green-wall and water/window
    3136! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
     
    33033308            DO  i = nxl, nxr
    33043309                DO  k = nzb_do, nzt_do
    3305 !                   print*, j,nys,nyn,i,nxl,nxr,k,nzb_do,nzt_do,local_pf(i,j,k),temp_pf(k,j,i)
    33063310                    local_pf(i,j,k) = temp_pf(k,j,i)
    33073311                ENDDO
Note: See TracChangeset for help on using the changeset viewer.