Ignore:
Timestamp:
Sep 30, 2019 8:40:37 AM (2 years ago)
Author:
pavelkrc
Message:

Merge branch resler into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r4226 r4245  
    2525! -----------------
    2626! $Id$
     27! Add reading and processing of building_surface_pars
     28!
     29! 4226 2019-09-10 17:03:24Z suehring
    2730! - Netcdf input routine for dimension length renamed
    2831! - Move offline-nesting-specific checks to nesting_offl_mod
     
    376379       REAL(wp)                              ::  fill1 = -9999.9_wp !< fill values for lod = 1
    377380       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d             !< 2d variable (lod = 1)
     381       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  oro_max            !< terraing height under particular buildings
    378382    END TYPE build_in
    379383
     
    416420       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  pars_xyz          !< respective parameters, level of detail = 2
    417421    END TYPE pars
     422!
     423!-- Data type for surface parameter lists
     424    TYPE pars_surf
     425       INTEGER(iwp)                                ::  np          !< total number of parameters
     426       INTEGER(iwp)                                ::  nsurf       !< number of local surfaces
     427       INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  index_ji    !< index for beginning and end of surfaces at (j,i)
     428       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  coords      !< (k,j,i,norm_z,norm_y,norm_x)
     429                                                                   !< k,j,i:                surface position
     430                                                                   !< norm_z,norm_y,norm_x: surface normal vector
     431
     432       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
     433
     434       REAL(wp)                              ::  fill = -9999.9_wp !< fill value
     435       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pars              !< respective parameters per surface
     436    END TYPE pars_surf
    418437!
    419438!-- Define type for global file attributes
     
    551570    TYPE(pars)  ::  water_pars_f               !< input variable for water parameters
    552571
     572    TYPE(pars_surf)  ::  building_surface_pars_f  !< input variable for building surface parameters
     573
    553574    TYPE(chem_emis_att_type)                             ::  chem_emis_att    !< Input Information of Chemistry Emission Data from netcdf 
    554575    TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:)  ::  chem_emis        !< Input Chemistry Emission Data from netcdf 
     
    683704!-- Public variables
    684705    PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f,    &
    685            building_id_f, building_pars_f, building_type_f,                    &
     706           building_id_f, building_pars_f, building_surface_pars_f,            &
     707           building_type_f,                                                    &
    686708           char_fill,                                                          &
    687709           char_lod,                                                           &
     
    19581980          building_pars_f%from_file = .FALSE.
    19591981       ENDIF
     1982!
     1983!--    Read building surface parameters
     1984       IF ( check_existence( var_names, 'building_surface_pars' ) )  THEN
     1985          building_surface_pars_f%from_file = .TRUE.
     1986          CALL get_attribute( id_surf, char_fill,                              &
     1987                              building_surface_pars_f%fill,                    &
     1988                              .FALSE., 'building_surface_pars' )
     1989!
     1990!--       Read building_surface_pars
     1991          CALL get_variable_surf( id_surf, 'building_surface_pars', &
     1992                                  building_surface_pars_f )
     1993       ELSE
     1994          building_surface_pars_f%from_file = .FALSE.
     1995       ENDIF
    19601996
    19611997!
     
    43024338             ENDIF
    43034339!
    4304 !--       Extrapolate if actual height is above the highest Inifor level
     4340!--       Extrapolate if actual height is above the highest Inifor level by the last value
    43054341          ELSE
    4306              var(k) = var_file(ku) + ( var_file(ku) - var_file(ku-1) ) /       &
    4307                                      ( z_file(ku)   - z_file(ku-1)   ) *       &
    4308                                      ( z_grid(k)    - z_file(ku)     )
    4309 
     4342             var(k) = var_file(ku)
    43104343          ENDIF
    43114344
     
    49654998#endif
    49664999    END SUBROUTINE get_variable_pr
     5000
     5001
     5002!------------------------------------------------------------------------------!
     5003! Description:
     5004! ------------
     5005!> Reads a per-surface pars variable from file. Because all surfaces are stored
     5006!> as flat 1-D array, each PE has to scan the data and find the surface indices
     5007!> belonging to its subdomain. During this scan, it also builds a necessary
     5008!> (j,i) index.
     5009!------------------------------------------------------------------------------!
     5010    SUBROUTINE get_variable_surf( id, variable_name, surf )
     5011#if defined( __netcdf )
     5012
     5013       USE pegrid
     5014
     5015       USE indices,                                            &
     5016           ONLY:  nxl, nxr, nys, nyn, nzb, nzt
     5017
     5018       USE control_parameters,                                 &
     5019           ONLY: dz, message_string
     5020
     5021       USE grid_variables,                                     &
     5022           ONLY: dx, dy
     5023       
     5024       USE basic_constants_and_equations_mod,                  &
     5025           ONLY: pi
     5026
     5027       USE arrays_3d,                                          &
     5028           ONLY:  zu, zw
     5029
     5030       IMPLICIT NONE
     5031
     5032       INTEGER, PARAMETER ::  nsurf_pars_read = 1024**2 !< read buffer size
     5033
     5034       CHARACTER(LEN=*)                          ::  variable_name !< variable name
     5035
     5036       INTEGER(iwp), DIMENSION(6)                ::  coords        !< integer coordinates of surface
     5037       INTEGER(iwp)                              ::  i, j
     5038       INTEGER(iwp)                              ::  isurf         !< netcdf surface index
     5039       INTEGER(iwp)                              ::  is            !< local surface index
     5040       INTEGER(iwp), INTENT(IN)                  ::  id            !< file id
     5041       INTEGER(iwp), DIMENSION(2)                ::  id_dim        !< dimension ids
     5042       INTEGER(iwp)                              ::  id_var        !< variable id
     5043       INTEGER(iwp)                              ::  id_zs         !< zs variable id
     5044       INTEGER(iwp)                              ::  id_ys         !< ys variable id
     5045       INTEGER(iwp)                              ::  id_xs         !< xs variable id
     5046       INTEGER(iwp)                              ::  id_zenith     !< zeith variable id
     5047       INTEGER(iwp)                              ::  id_azimuth    !< azimuth variable id
     5048       INTEGER(iwp)                              ::  is0, isc      !< read surface start and count
     5049       INTEGER(iwp)                              ::  nsurf         !< total number of surfaces in file
     5050       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nsurf_ji      !< numbers of surfaces by coords
     5051
     5052       TYPE(pars_surf)                           ::  surf          !< parameters variable to be loaded
     5053       REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  pars_read     !< read buffer
     5054       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  zs, ys, xs    !< read buffer for zs(s), ys, xs
     5055       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  zenith        !< read buffer for zenith(s)
     5056       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  azimuth       !< read buffer for azimuth(s)
     5057       REAL(wp)                                  ::  oro_max_l     !< maximum terrain height under building
     5058
     5059!
     5060!--    First, inquire variable ID
     5061       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
     5062       nc_stat = NF90_INQ_VARID( id, 'zs',                  id_zs )
     5063       nc_stat = NF90_INQ_VARID( id, 'ys',                  id_ys )
     5064       nc_stat = NF90_INQ_VARID( id, 'xs',                  id_xs )
     5065       nc_stat = NF90_INQ_VARID( id, 'zenith',              id_zenith )
     5066       nc_stat = NF90_INQ_VARID( id, 'azimuth',             id_azimuth )
     5067!
     5068!--    Inquire dimension sizes
     5069       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
     5070       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = nsurf )
     5071       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(2), LEN = surf%np )
     5072
     5073       ALLOCATE ( pars_read( nsurf_pars_read, surf%np ),        &
     5074                  zs(nsurf_pars_read), ys(nsurf_pars_read),     &
     5075                  xs(nsurf_pars_read), zenith(nsurf_pars_read), &
     5076                  azimuth(nsurf_pars_read),                     &
     5077                  nsurf_ji(nys:nyn, nxl:nxr) )
     5078
     5079       nsurf_ji(:,:) = 0
     5080!
     5081!--    Scan surface coordinates, count local
     5082       is0 = 1
     5083       DO
     5084          isc = MIN(nsurf_pars_read, nsurf - is0 + 1)
     5085          IF ( isc <= 0 )  EXIT
     5086
     5087          nc_stat = NF90_GET_VAR( id, id_ys, ys,     &
     5088                                  start = (/ is0 /), &
     5089                                  count = (/ isc /) )
     5090          nc_stat = NF90_GET_VAR( id, id_xs, xs,     &
     5091                                  start = (/ is0 /), &
     5092                                  count = (/ isc /) )
     5093          nc_stat = NF90_GET_VAR( id, id_zenith, zenith,      &
     5094                                  start = (/ is0 /), &
     5095                                  count = (/ isc /) )
     5096          nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth,    &
     5097                                  start = (/ is0 /), &
     5098                                  count = (/ isc /) )
     5099          CALL handle_error( 'get_variable_surf', 682, 'azimuth' )
     5100         
     5101          DO  isurf = 1, isc
     5102!
     5103!--          Parse coordinates, detect if belongs to subdomain
     5104             coords = transform_coords( xs(isurf), ys(isurf),         &
     5105                                        zenith(isurf), azimuth(isurf) )
     5106             IF ( coords(2) < nys  .OR.  coords(2) > nyn  .OR.  &
     5107                  coords(3) < nxl  .OR.  coords(3) > nxr )  CYCLE
     5108
     5109             nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1
     5110          ENDDO
     5111
     5112          is0 = is0 + isc
     5113       ENDDO
     5114!
     5115!--    Populate reverse index from surface counts
     5116       ALLOCATE ( surf%index_ji( 2, nys:nyn, nxl:nxr ) )
     5117
     5118       isurf = 1
     5119       DO  j = nys, nyn
     5120          DO  i = nxl, nxr
     5121             surf%index_ji(:,j,i) = (/ isurf, isurf + nsurf_ji(j,i) - 1 /)
     5122             isurf = isurf + nsurf_ji(j,i)
     5123          ENDDO
     5124       ENDDO
     5125
     5126       surf%nsurf = isurf - 1
     5127       ALLOCATE( surf%pars( 0:surf%np-1, surf%nsurf ), &
     5128                 surf%coords( 6, surf%nsurf ) )
     5129!
     5130!--    Scan surfaces again, saving pars into allocated structures
     5131       nsurf_ji(:,:) = 0
     5132       is0 = 1
     5133       DO
     5134          isc = MIN(nsurf_pars_read, nsurf - is0 + 1)
     5135          IF ( isc <= 0 )  EXIT
     5136
     5137          nc_stat = NF90_GET_VAR( id, id_var, pars_read(1:isc, 1:surf%np), &
     5138                                  start = (/ is0, 1       /),              &
     5139                                  count = (/ isc, surf%np /) )
     5140          CALL handle_error( 'get_variable_surf', 683, variable_name )
     5141
     5142          nc_stat = NF90_GET_VAR( id, id_zs, zs,                           &
     5143                                  start = (/ is0 /),                       &
     5144                                  count = (/ isc /) )
     5145          nc_stat = NF90_GET_VAR( id, id_ys, ys,                           &
     5146                                  start = (/ is0 /),                       &
     5147                                  count = (/ isc /) )
     5148          nc_stat = NF90_GET_VAR( id, id_xs, xs,                           &
     5149                                  start = (/ is0 /),                       &
     5150                                  count = (/ isc /) )
     5151          nc_stat = NF90_GET_VAR( id, id_zenith, zenith,                   &
     5152                                  start = (/ is0 /),                       &
     5153                                  count = (/ isc /) )
     5154          nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth,                 &
     5155                                  start = (/ is0 /),                       &
     5156                                  count = (/ isc /) )
     5157         
     5158          DO  isurf = 1, isc
     5159!
     5160!--          Parse coordinates, detect if belongs to subdomain
     5161             coords = transform_coords( xs(isurf), ys(isurf),         &
     5162                                        zenith(isurf), azimuth(isurf) )
     5163             IF ( coords(2) < nys  .OR.  coords(2) > nyn  .OR.  &
     5164                  coords(3) < nxl  .OR.  coords(3) > nxr )  CYCLE
     5165!
     5166!--          Determine maximum terrain under building (base z-coordinate). Using
     5167!--          normal vector to locate building inner coordinates.
     5168             oro_max_l = buildings_f%oro_max(coords(2)-coords(5), coords(3)-coords(6))
     5169             IF  ( oro_max_l == buildings_f%fill1 )  THEN
     5170                WRITE( message_string, * ) 'Found building surface on '   // &
     5171                   'non-building coordinates (xs, ys, zenith, azimuth): ',   &
     5172                   xs(isurf), ys(isurf), zenith(isurf), azimuth(isurf)
     5173                CALL message( 'get_variable_surf', 'PA0684', 2, 2, 0, 6, 0 )
     5174             ENDIF
     5175!
     5176!--          Urban layer has no stretching, therefore using dz(1) instead of linear
     5177!--          searching through zu/zw
     5178             coords(1) = NINT((zs(isurf) + oro_max_l) / dz(1) +     &
     5179                              0.5_wp + 0.5_wp * coords(4), KIND=iwp)
     5180!
     5181!--          Save surface entry
     5182             is = surf%index_ji(1, coords(2), coords(3)) + nsurf_ji(coords(2), coords(3))
     5183             surf%pars(:,is) = pars_read(isurf,:)
     5184             surf%coords(:,is) = coords(:)
     5185
     5186             nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1
     5187          ENDDO
     5188
     5189          is0 = is0 + isc
     5190       ENDDO
     5191
     5192       DEALLOCATE( pars_read, zs, ys, xs, zenith, azimuth, nsurf_ji )
     5193
     5194    CONTAINS
     5195
     5196       PURE FUNCTION transform_coords( x, y, zenith, azimuth )
     5197
     5198          REAL(wp), INTENT(in)       ::  x, y    !< surface centre coordinates in metres from origin
     5199          REAL(wp), INTENT(in)       ::  zenith  !< surface normal zenith angle in degrees
     5200          REAL(wp), INTENT(in)       ::  azimuth !< surface normal azimuth angle in degrees
     5201
     5202          INTEGER(iwp), DIMENSION(6) ::  transform_coords !< (k,j,i,norm_z,norm_y,norm_x)
     5203
     5204          transform_coords(4) = NINT(COS(zenith*pi/180._wp), KIND=iwp)
     5205          IF ( transform_coords(4) == 0 )  THEN
     5206             transform_coords(5) = NINT(COS(azimuth*pi/180._wp), KIND=iwp)
     5207             transform_coords(6) = NINT(SIN(azimuth*pi/180._wp), KIND=iwp)
     5208          ELSE
     5209             transform_coords(5) = 0._wp
     5210             transform_coords(6) = 0._wp
     5211          ENDIF
     5212
     5213          transform_coords(1) = -999._wp ! not calculated here
     5214          transform_coords(2) = NINT(y/dy - 0.5_wp + 0.5_wp*transform_coords(5), KIND=iwp)
     5215          transform_coords(3) = NINT(x/dx - 0.5_wp + 0.5_wp*transform_coords(6), KIND=iwp)
     5216
     5217       END FUNCTION transform_coords
     5218
     5219#endif
     5220    END SUBROUTINE get_variable_surf
    49675221
    49685222
Note: See TracChangeset for help on using the changeset viewer.