Changeset 3182 for palm/trunk/UTIL/inifor/src/grid.f90
- Timestamp:
- Jul 27, 2018 1:36:03 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/UTIL/inifor/src/grid.f90 (modified) (93 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/grid.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! 23 ! Introduced new PALM grid stretching 24 ! Updated variable names and metadata for PIDS v1.9 compatibility 25 ! Better compatibility with older Intel compilers: 26 ! - avoiding implicit array allocation with new get_netcdf_variable() 27 ! subroutine instead of function 28 ! Improved command line interface: 29 ! - Produce forcing variables when '-mode profile' is being used 30 ! - Renamend initial-condition mode variable 'mode' to 'ic_mode' 31 ! - Improved handling of the start date string 32 ! Removed unnecessary variables and routines 33 ! 24 34 ! 25 35 ! Former revisions: … … 47 57 USE defs, & 48 58 ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp, & 49 SNAME, LNAME, PATH, FORCING_ FREQ, WATER_ID, FILL_ITERATIONS, &59 SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, & 50 60 BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L 51 61 USE io, & 52 62 ONLY: get_netcdf_variable, get_netcdf_attribute, & 53 parse_command_line_arguments 63 parse_command_line_arguments, validate_config 54 64 USE netcdf, & 55 65 ONLY: NF90_MAX_NAME, NF90_MAX_VAR_DIMS … … 67 77 SAVE 68 78 69 REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad] 70 REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad] 71 REAL(dp) :: lambda_n = 0.0_dp !< longitude of rotaded pole of COSMO-DE grid [rad] 72 REAL(dp) :: phi_c = 0.0_dp !< rotated-grid latitude of the center of the PALM domain [rad] 73 REAL(dp) :: lambda_c = 0.0_dp !< rotated-grid longitude of the centre of the PALM domain [rad] 74 REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] 75 REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] 76 REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug 77 REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m] 78 REAL(dp) :: dy = 0.0_dp !< PALM-4U grid spacing in y direction [m] 79 REAL(dp) :: dz = 0.0_dp !< PALM-4U grid spacing in z direction [m] 80 REAL(dp) :: dxi = 0.0_dp !< inverse PALM-4U grid spacing in x direction [m^-1] 81 REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1] 82 REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1] 83 REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m] 84 REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m] 85 REAL(dp) :: lz = 0.0_dp !< PALM-4U domain size in z direction [m] 86 REAL(dp) :: ug = 0.0_dp !< geostrophic wind in x direction [m/s] 87 REAL(dp) :: vg = 0.0_dp !< geostrophic wind in y direction [m/s] 88 REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa] 89 REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m] 90 REAL(dp) :: y0 = 0.0_dp !< y coordinate of PALM-4U Earth tangent [m] 91 REAL(dp) :: z0 = 0.0_dp !< Elevation of the PALM-4U domain above sea level [m] 92 REAL(dp) :: lonmin = 0.0_dp !< Minimunm longitude of COSMO-DE's rotated-pole grid 93 REAL(dp) :: lonmax = 0.0_dp !< Maximum longitude of COSMO-DE's rotated-pole grid 94 REAL(dp) :: latmin = 0.0_dp !< Minimunm latitude of COSMO-DE's rotated-pole grid 95 REAL(dp) :: latmax = 0.0_dp !< Maximum latitude of COSMO-DE's rotated-pole grid 96 REAL(dp) :: latitude = 0.0_dp !< geograpohical latitude of the PALM-4U origin, from inipar namelist [deg] 97 REAL(dp) :: longitude = 0.0_dp !< geograpohical longitude of the PALM-4U origin, from inipar namelist [deg] 98 REAL(dp) :: origin_lat= 0.0_dp !< geograpohical latitude of the PALM-4U origin, from static driver netCDF file [deg] 99 REAL(dp) :: origin_lon= 0.0_dp !< geograpohical longitude of the PALM-4U origin, from static driver netCDF file [deg] 100 REAL(dp) :: end_time = 0.0_dp !< PALM-4U simulation time [s] 79 REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad] 80 REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad] 81 REAL(dp) :: lambda_n = 0.0_dp !< longitude of rotaded pole of COSMO-DE grid [rad] 82 REAL(dp) :: phi_c = 0.0_dp !< rotated-grid latitude of the center of the PALM domain [rad] 83 REAL(dp) :: lambda_c = 0.0_dp !< rotated-grid longitude of the centre of the PALM domain [rad] 84 REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] 85 REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] 86 REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug 87 REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m] 88 REAL(dp) :: dy = 0.0_dp !< PALM-4U grid spacing in y direction [m] 89 REAL(dp) :: dz(10) = -1.0_dp !< PALM-4U grid spacing in z direction [m] 90 REAL(dp) :: dz_max = 1000.0_dp !< maximum vertical grid spacing [m] 91 REAL(dp) :: dz_stretch_factor = 1.08_dp !< factor for vertical grid stretching [m] 92 REAL(dp) :: dz_stretch_level = -9999999.9_dp!< height above which the vertical grid will be stretched [m] 93 REAL(dp) :: dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter 94 REAL(dp) :: dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter 95 REAL(dp) :: dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter 96 REAL(dp) :: dxi = 0.0_dp !< inverse PALM-4U grid spacing in x direction [m^-1] 97 REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1] 98 REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1] 99 REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m] 100 REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m] 101 REAL(dp) :: ug = 0.0_dp !< geostrophic wind in x direction [m/s] 102 REAL(dp) :: vg = 0.0_dp !< geostrophic wind in y direction [m/s] 103 REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa] 104 REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m] 105 REAL(dp) :: y0 = 0.0_dp !< y coordinate of PALM-4U Earth tangent [m] 106 REAL(dp) :: z0 = 0.0_dp !< Elevation of the PALM-4U domain above sea level [m] 107 REAL(dp) :: z_top = 0.0_dp !< height of the scalar top boundary [m] 108 REAL(dp) :: zw_top = 0.0_dp !< height of the vertical velocity top boundary [m] 109 REAL(dp) :: lonmin = 0.0_dp !< Minimunm longitude of COSMO-DE's rotated-pole grid 110 REAL(dp) :: lonmax = 0.0_dp !< Maximum longitude of COSMO-DE's rotated-pole grid 111 REAL(dp) :: latmin = 0.0_dp !< Minimunm latitude of COSMO-DE's rotated-pole grid 112 REAL(dp) :: latmax = 0.0_dp !< Maximum latitude of COSMO-DE's rotated-pole grid 113 REAL(dp) :: latitude = 0.0_dp !< geographical latitude of the PALM-4U origin, from inipar namelist [deg] 114 REAL(dp) :: longitude = 0.0_dp !< geographical longitude of the PALM-4U origin, from inipar namelist [deg] 115 REAL(dp) :: origin_lat = 0.0_dp !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg] 116 REAL(dp) :: origin_lon = 0.0_dp !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg] 117 REAL(dp) :: rotation_angle = 0.0_dp !< clockwise angle the PALM-4U north is rotated away from geographical north [deg] 118 REAL(dp) :: end_time = 0.0_dp !< PALM-4U simulation time [s] 101 119 102 120 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file … … 106 124 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlon !< longitudes of COSMO-DE's rotated-pole grid 107 125 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlat !< latitudes of COSMO-DE's rotated-pole grid 108 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: time 126 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: time !< output times 127 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: x !< base palm grid x coordinate vector pointed to by grid_definitions 128 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: xu !< base palm grid xu coordinate vector pointed to by grid_definitions 129 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: y !< base palm grid y coordinate vector pointed to by grid_definitions 130 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: yv !< base palm grid yv coordinate vector pointed to by grid_definitions 131 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z_column !< base palm grid z coordinate vector including the top boundary coordinate (entire column) 132 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw_column !< base palm grid zw coordinate vector including the top boundary coordinate (entire column) 133 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z !< base palm grid z coordinate vector pointed to by grid_definitions 134 REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw !< base palm grid zw coordinate vector pointed to by grid_definitions 109 135 110 136 INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: soiltyp !< COSMO-DE soil type map 137 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 138 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 111 139 INTEGER :: i !< indexing variable 112 INTEGER :: imin, imax, jmin,jmax !< index ranges for profile averaging140 INTEGER :: average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging 113 141 INTEGER :: k !< indexing variable 114 142 INTEGER :: nt !< number of output time steps … … 130 158 LOGICAL :: init_variables_required 131 159 LOGICAL :: boundary_variables_required 160 LOGICAL :: ls_forcing_variables_required 161 LOGICAL :: profile_grids_required 132 162 133 163 INTEGER :: n_invar = 0 !< number of variables in the input variable table … … 193 223 TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable 194 224 195 NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude 225 NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude, & 226 dz_max, dz_stretch_factor, dz_stretch_level, & !< old grid stretching parameters 227 dz_stretch_level_start, dz_stretch_level_end !< new grid stretching parameters 196 228 NAMELIST /d3par/ end_time 197 229 198 230 CHARACTER(LEN=LNAME) :: nc_source_text = '' !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...' 199 CHARACTER(LEN=DATE) :: start_date = '' !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation200 CHARACTER(LEN=PATH) :: hhl_file = '' !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces)201 CHARACTER(LEN=PATH) :: namelist_file = '' !< Path to the PALM-4U namelist file202 CHARACTER(LEN=PATH) :: static_driver_file = '' !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)203 CHARACTER(LEN=PATH) :: soiltyp_file = '' !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types)204 CHARACTER(LEN=PATH) :: input_path = '' !< Path to the input data file directory205 231 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: flow_files 206 232 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_moisture_files … … 212 238 CHARACTER(LEN=SNAME) :: radiation_suffix !< Suffix of radiation input files, e.g. 'radiation' 213 239 CHARACTER(LEN=SNAME) :: soilmoisture_suffix !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture' 214 CHARACTER(LEN=SNAME) :: mode !< INIFOR's initialization mode, 'profile' or 'volume'215 240 216 241 TYPE(nc_file) :: output_file 242 243 TYPE(inifor_config) :: cfg !< Container of INIFOR configuration 217 244 218 245 CONTAINS … … 224 251 ! Section 1: Define default parameters 225 252 !------------------------------------------------------------------------------ 226 start_date = '2013072100'253 cfg % start_date = '2013072100' 227 254 end_hour = 2 228 255 start_hour_soil = -2 … … 249 276 origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation 250 277 origin_lon = 13.082744_dp * TO_RADIANS 251 z0 = 35.0_dp278 cfg % z0 = 35.0_dp 252 279 253 280 ! Default atmospheric parameters 254 281 ug = 0.0_dp 255 282 vg = 0.0_dp 256 p0 = P_SL283 cfg % p0 = P_SL 257 284 258 285 ! Parameters for file names … … 262 289 start_hour_soilmoisture = start_hour_flow - 2 263 290 end_hour_soilmoisture = start_hour_flow 264 step_hour = 1 291 step_hour = FORCING_STEP 292 265 293 input_prefix = 'laf' ! 'laf' for COSMO-DE analyses 294 cfg % flow_prefix = input_prefix 295 cfg % soil_prefix = input_prefix 296 cfg % radiation_prefix = input_prefix 297 cfg % soilmoisture_prefix = input_prefix 298 266 299 flow_suffix = '-flow' 267 300 soil_suffix = '-soil' … … 273 306 !------------------------------------------------------------------------------ 274 307 275 ! Set default paths 276 input_path = './' 277 hhl_file = '' 278 soiltyp_file = '' 279 namelist_file = './namelist' 280 output_file % name = './palm-4u-input.nc' 308 ! Set default paths and modes 309 cfg % input_path = './' 310 cfg % hhl_file = '' 311 cfg % soiltyp_file = '' 312 cfg % namelist_file = './namelist' 313 cfg % static_driver_file = '' 314 cfg % output_file = './palm-4u-input.nc' 315 cfg % ic_mode = 'volume' 316 cfg % bc_mode = 'real' 281 317 282 318 ! Update default file names and domain centre 283 CALL parse_command_line_arguments( start_date, hhl_file, soiltyp_file, & 284 static_driver_file, input_path, output_file % name, & 285 namelist_file, ug, vg, p0, z0, mode ) 286 287 init_variables_required = .TRUE. 288 boundary_variables_required = (TRIM(mode) .NE. 'profile') 289 290 CALL normalize_path(input_path) 291 IF (TRIM(hhl_file) == '') hhl_file = TRIM(input_path) // 'hhl.nc' 292 IF (TRIM(soiltyp_file) == '') soiltyp_file = TRIM(input_path) // 'soil.nc' 293 294 CALL report('setup_parameters', " data path: " // TRIM(input_path)) 295 CALL report('setup_parameters', " hhl file: " // TRIM(hhl_file)) 296 CALL report('setup_parameters', " soiltyp file: " // TRIM(soiltyp_file)) 297 CALL report('setup_parameters', " namelist file: " // TRIM(namelist_file)) 298 CALL report('setup_parameters', "output data file: " // TRIM(output_file % name)) 319 CALL parse_command_line_arguments( cfg ) 320 321 output_file % name = cfg % output_file 322 z0 = cfg % z0 323 p0 = cfg % p0 324 325 init_variables_required = .TRUE. 326 boundary_variables_required = TRIM( cfg % bc_mode ) == 'real' 327 ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal' 328 329 IF ( ls_forcing_variables_required ) THEN 330 message = "Averaging of large-scale forcing profiles " // & 331 "has not been implemented, yet." 332 CALL abort('setup_parameters', message) 333 END IF 334 335 CALL normalize_path(cfg % input_path) 336 IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' 337 IF (TRIM(cfg % soiltyp_file) == '') cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc' 338 339 CALL validate_config( cfg ) 340 341 CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode)) 342 CALL report('setup_parameters', " forcing mode: " // TRIM(cfg % bc_mode)) 343 CALL report('setup_parameters', " data path: " // TRIM(cfg % input_path)) 344 CALL report('setup_parameters', " hhl file: " // TRIM(cfg % hhl_file)) 345 CALL report('setup_parameters', " soiltyp file: " // TRIM(cfg % soiltyp_file)) 346 CALL report('setup_parameters', " namelist file: " // TRIM(cfg % namelist_file)) 347 CALL report('setup_parameters', " output data file: " // TRIM(output_file % name)) 299 348 300 349 CALL run_control('time', 'init') 301 350 ! Read in namelist parameters 302 OPEN(10, FILE= namelist_file)351 OPEN(10, FILE=cfg % namelist_file) 303 352 READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz 304 353 READ(10, NML=d3par) ! end_time … … 306 355 CALL run_control('time', 'read') 307 356 308 end_hour = CEILING( end_time / FORCING_FREQ)357 end_hour = CEILING( end_time / 3600.0 * step_hour ) 309 358 310 359 ! Generate input file lists 311 CALL input_file_list(start_date, start_hour_flow, end_hour, step_hour, &312 input_path, input_prefix, flow_suffix, flow_files)313 CALL input_file_list(start_date, start_hour_soil, end_hour, step_hour, &314 input_path, input_prefix, soil_suffix, soil_files)315 CALL input_file_list(start_date, start_hour_radiation, end_hour, step_hour, &316 input_path, input_prefix, radiation_suffix, radiation_files)317 CALL input_file_list(start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &318 input_path, input_prefix, soilmoisture_suffix, soil_moisture_files)360 CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour, & 361 cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files) 362 CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour, & 363 cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files) 364 CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, & 365 cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files) 366 CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 367 cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) 319 368 320 369 ! … … 322 371 ! Section 3: Check for consistency 323 372 !------------------------------------------------------------------------------ 324 IF (dx*dy*dz .EQ. 0.0_dp) THEN 325 message = "Grid cells have zero volume. Grid spacings are probably"//& 326 " set incorrectly in namelist file '" // TRIM(namelist_file) // "'." 327 CALL abort('setup_parameters', message) 328 END IF 373 329 374 ! 330 375 !------------------------------------------------------------------------------ … … 339 384 ! Read COSMO-DE soil type map 340 385 cosmo_var % name = 'SOILTYP' 341 soiltyp = NINT(get_netcdf_variable(soiltyp_file, cosmo_var), hp)386 CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp) 342 387 343 388 message = 'Reading PALM-4U origin from' 344 IF (TRIM( static_driver_file) .NE. '') THEN345 346 origin_lon = get_netcdf_attribute( static_driver_file, 'origin_lon')347 origin_lat = get_netcdf_attribute( static_driver_file, 'origin_lat')389 IF (TRIM(cfg % static_driver_file) .NE. '') THEN 390 391 origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon') 392 origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat') 348 393 349 394 message = TRIM(message) // " static driver file '" & 350 // TRIM( static_driver_file) // "'"395 // TRIM(cfg % static_driver_file) // "'" 351 396 352 397 … … 357 402 358 403 message = TRIM(message) // " namlist file '" & 359 // TRIM( namelist_file) // "'"404 // TRIM(cfg % namelist_file) // "'" 360 405 361 406 END IF … … 368 413 ! Read in COSMO-DE heights of half layers (vertical cell faces) 369 414 cosmo_var % name = 'HHL' 370 hhl = get_netcdf_variable(hhl_file, cosmo_var)415 CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl) 371 416 CALL run_control('time', 'read') 372 417 … … 392 437 lx = (nx+1) * dx 393 438 ly = (ny+1) * dy 394 lz = (nz+1) * dz395 439 396 440 ! PALM-4U point of Earth tangency … … 399 443 400 444 ! time vector 401 nt = CEILING(end_time / FORCING_FREQ) + 1445 nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1 402 446 ALLOCATE( time(nt) ) 403 447 CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time) … … 463 507 SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...)) 464 508 CHARACTER :: interp_mode 465 509 466 510 !------------------------------------------------------------------------------ 467 ! Section 1: Define model and initialization grids 511 ! Section 0: Define base PALM-4U grid coordinate vectors 512 !------------------------------------------------------------------------------ 513 ! palm x y z, we allocate the column to nz+1 in order to include the top 514 ! scalar boundary. The interpolation grids will be associated with 515 ! a shorter column that omits the top element. 516 517 ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) ) 518 CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x) 519 CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y) 520 CALL stretched_z(z_column, dz, dz_max=dz_max, & 521 dz_stretch_factor=dz_stretch_factor, & 522 dz_stretch_level=dz_stretch_level, & 523 dz_stretch_level_start=dz_stretch_level_start, & 524 dz_stretch_level_end=dz_stretch_level_end, & 525 dz_stretch_factor_array=dz_stretch_factor_array) 526 z(1:nz) = z_column(1:nz) 527 z_top = z_column(nz+1) 528 529 ! palm xu yv zw, compared to the scalar grid, velocity coordinates 530 ! contain one element less. 531 ALLOCATE( xu(1:nx), yv(1:ny), zw(1:nz-1), zw_column(1:nz)) 532 CALL linspace(dx, lx - dx, xu) 533 CALL linspace(dy, ly - dy, yv) 534 CALL midpoints(z_column, zw_column) 535 zw(1:nz-1) = zw_column(1:nz-1) 536 zw_top = zw_column(nz) 537 538 539 !------------------------------------------------------------------------------ 540 ! Section 1: Define initialization and boundary grids 468 541 !------------------------------------------------------------------------------ 469 542 CALL init_grid_definition('palm', grid=palm_grid, & 470 543 xmin=0.0_dp, xmax=lx, & 471 544 ymin=0.0_dp, ymax=ly, & 472 zmin=0.0_dp, zmax=lz,x0=x0, y0=y0, z0=z0, &473 nx=nx, ny=ny, nz=nz, mode=mode)545 x0=x0, y0=y0, z0=z0, & 546 nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode) 474 547 475 548 ! Subtracting 1 because arrays will be allocated with nlon + 1 elements. … … 477 550 xmin=lonmin, xmax=lonmax, & 478 551 ymin=latmin, ymax=latmax, & 479 zmin=0.0_dp, zmax=51.0_dp,x0=x0, y0=y0, z0=0.0_dp, &552 x0=x0, y0=y0, z0=0.0_dp, & 480 553 nx=nlon-1, ny=nlat-1, nz=nlev-1) 481 554 … … 487 560 xmin=0.0_dp, xmax=lx, & 488 561 ymin=0.0_dp, ymax=ly, & 489 zmin=0.0_dp, zmax=lz,x0=x0, y0=y0, z0=z0, &562 x0=x0, y0=y0, z0=z0, & 490 563 nx=nx, ny=ny, nz=nlev-2) 491 564 … … 493 566 xmin = dx, xmax = lx - dx, & 494 567 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 495 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &496 568 x0=x0, y0=y0, z0 = z0, & 497 569 nx = nx-1, ny = ny, nz = nz, & 498 dx = dx, dy = dy, dz = dz, mode=mode)570 z=z, ic_mode=cfg % ic_mode) 499 571 500 572 CALL init_grid_definition('boundary', grid=v_initial_grid, & 501 573 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 502 574 ymin = dy, ymax = ly - dy, & 503 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &504 575 x0=x0, y0=y0, z0 = z0, & 505 576 nx = nx, ny = ny-1, nz = nz, & 506 dx = dx, dy = dy, dz = dz, mode=mode)577 z=z, ic_mode=cfg % ic_mode) 507 578 508 579 CALL init_grid_definition('boundary', grid=w_initial_grid, & 509 580 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 510 581 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 511 zmin = dz, zmax = lz - dz, &512 582 x0=x0, y0=y0, z0 = z0, & 513 583 nx = nx, ny = ny, nz = nz-1, & 514 dx = dx, dy = dy, dz = dz, mode=mode)584 z=zw, ic_mode=cfg % ic_mode) 515 585 516 586 CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate, & 517 587 xmin = dx, xmax = lx - dx, & 518 588 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 519 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &520 589 x0=x0, y0=y0, z0 = z0, & 521 nx = nx-1, ny = ny, nz = nlev - 2, & 522 dx = dx, dy = dy, dz = dz) 590 nx = nx-1, ny = ny, nz = nlev - 2) 523 591 524 592 CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate, & 525 593 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 526 594 ymin = dy, ymax = ly - dy, & 527 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &528 595 x0=x0, y0=y0, z0 = z0, & 529 nx = nx, ny = ny-1, nz = nlev - 2, & 530 dx = dx, dy = dy, dz = dz) 596 nx = nx, ny = ny-1, nz = nlev - 2) 531 597 532 598 CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate, & 533 599 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 534 600 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 535 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &536 601 x0=x0, y0=y0, z0 = z0, & 537 nx = nx, ny = ny, nz = nlev - 2, & 538 dx = dx, dy = dy, dz = dz) 602 nx = nx, ny = ny, nz = nlev - 2) 539 603 540 604 IF (boundary_variables_required) THEN … … 546 610 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 547 611 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 548 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 549 x0=x0, y0=y0, z0 = z0, & 550 nx = 0, ny = ny, nz = nz, & 551 dx = dx, dy = dy, dz = dz) 612 x0=x0, y0=y0, z0 = z0, & 613 nx = 0, ny = ny, nz = nz, z=z) 552 614 553 615 CALL init_grid_definition('boundary', grid=scalars_west_grid, & 554 616 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 555 617 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 556 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 557 x0=x0, y0=y0, z0 = z0, & 558 nx = 0, ny = ny, nz = nz, & 559 dx = dx, dy = dy, dz = dz) 618 x0=x0, y0=y0, z0 = z0, & 619 nx = 0, ny = ny, nz = nz, z=z) 560 620 561 621 CALL init_grid_definition('boundary', grid=scalars_north_grid, & 562 622 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 563 623 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 564 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 565 x0=x0, y0=y0, z0 = z0, & 566 nx = nx, ny = 0, nz = nz, & 567 dx = dx, dy = dy, dz = dz) 624 x0=x0, y0=y0, z0 = z0, & 625 nx = nx, ny = 0, nz = nz, z=z) 568 626 569 627 CALL init_grid_definition('boundary', grid=scalars_south_grid, & 570 628 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 571 629 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 572 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 573 x0=x0, y0=y0, z0 = z0, & 574 nx = nx, ny = 0, nz = nz, & 575 dx = dx, dy = dy, dz = dz) 630 x0=x0, y0=y0, z0 = z0, & 631 nx = nx, ny = 0, nz = nz, z=z) 576 632 577 633 CALL init_grid_definition('boundary', grid=scalars_top_grid, & 578 634 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 579 635 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 580 zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz, & 581 x0=x0, y0=y0, z0 = z0, & 582 nx = nx, ny = ny, nz = 0, & 583 dx = dx, dy = dy, dz = dz) 636 x0=x0, y0=y0, z0 = z0, & 637 nx = nx, ny = ny, nz = 1, z=(/z_top/)) 584 638 585 639 CALL init_grid_definition('boundary', grid=u_east_grid, & 586 640 xmin = lx, xmax = lx, & 587 641 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 588 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 589 x0=x0, y0=y0, z0 = z0, & 590 nx = 0, ny = ny, nz = nz, & 591 dx = dx, dy = dy, dz = dz) 642 x0=x0, y0=y0, z0 = z0, & 643 nx = 0, ny = ny, nz = nz, z=z) 592 644 593 645 CALL init_grid_definition('boundary', grid=u_west_grid, & 594 646 xmin = 0.0_dp, xmax = 0.0_dp, & 595 647 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 596 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 597 x0=x0, y0=y0, z0 = z0, & 598 nx = 0, ny = ny, nz = nz, & 599 dx = dx, dy = dy, dz = dz) 648 x0=x0, y0=y0, z0 = z0, & 649 nx = 0, ny = ny, nz = nz, z=z) 600 650 601 651 CALL init_grid_definition('boundary', grid=u_north_grid, & 602 652 xmin = dx, xmax = lx - dx, & 603 653 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 604 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 605 x0=x0, y0=y0, z0 = z0, & 606 nx = nx-1, ny = 0, nz = nz, & 607 dx = dx, dy = dy, dz = dz) 608 654 x0=x0, y0=y0, z0 = z0, & 655 nx = nx-1, ny = 0, nz = nz, z=z) 656 609 657 CALL init_grid_definition('boundary', grid=u_south_grid, & 610 658 xmin = dx, xmax = lx - dx, & 611 659 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 612 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 613 x0=x0, y0=y0, z0 = z0, & 614 nx = nx-1, ny = 0, nz = nz, & 615 dx = dx, dy = dy, dz = dz) 660 x0=x0, y0=y0, z0 = z0, & 661 nx = nx-1, ny = 0, nz = nz, z=z) 616 662 617 663 CALL init_grid_definition('boundary', grid=u_top_grid, & 618 664 xmin = dx, xmax = lx - dx, & 619 665 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 620 zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz, & 621 x0=x0, y0=y0, z0 = z0, & 622 nx = nx-1, ny = ny, nz = 0, & 623 dx = dx, dy = dy, dz = dz) 666 x0=x0, y0=y0, z0 = z0, & 667 nx = nx-1, ny = ny, nz = 1, z=(/z_top/)) 624 668 625 669 CALL init_grid_definition('boundary', grid=v_east_grid, & 626 670 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 627 671 ymin = dy, ymax = ly - dy, & 628 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 629 x0=x0, y0=y0, z0 = z0, & 630 nx = 0, ny = ny-1, nz = nz, & 631 dx = dx, dy = dy, dz = dz) 672 x0=x0, y0=y0, z0 = z0, & 673 nx = 0, ny = ny-1, nz = nz, z=z) 632 674 633 675 CALL init_grid_definition('boundary', grid=v_west_grid, & 634 676 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 635 677 ymin = dy, ymax = ly - dy, & 636 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 637 x0=x0, y0=y0, z0 = z0, & 638 nx = 0, ny = ny-1, nz = nz, & 639 dx = dx, dy = dy, dz = dz) 678 x0=x0, y0=y0, z0 = z0, & 679 nx = 0, ny = ny-1, nz = nz, z=z) 640 680 641 681 CALL init_grid_definition('boundary', grid=v_north_grid, & 642 682 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 643 683 ymin = ly, ymax = ly, & 644 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 645 x0=x0, y0=y0, z0 = z0, & 646 nx = nx, ny = 0, nz = nz, & 647 dx = dx, dy = dy, dz = dz) 684 x0=x0, y0=y0, z0 = z0, & 685 nx = nx, ny = 0, nz = nz, z=z) 648 686 649 687 CALL init_grid_definition('boundary', grid=v_south_grid, & 650 688 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 651 689 ymin = 0.0_dp, ymax = 0.0_dp, & 652 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, & 653 x0=x0, y0=y0, z0 = z0, & 654 nx = nx, ny = 0, nz = nz, & 655 dx = dx, dy = dy, dz = dz) 690 x0=x0, y0=y0, z0 = z0, & 691 nx = nx, ny = 0, nz = nz, z=z) 656 692 657 693 CALL init_grid_definition('boundary', grid=v_top_grid, & 658 694 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 659 695 ymin = dy, ymax = ly - dy, & 660 zmin = lz + 0.5_dp * dz, zmax = lz + 0.5_dp * dz, & 661 x0=x0, y0=y0, z0 = z0, & 662 nx = nx, ny = ny-1, nz = 0, & 663 dx = dx, dy = dy, dz = dz) 696 x0=x0, y0=y0, z0 = z0, & 697 nx = nx, ny = ny-1, nz = 1, z=(/z_top/)) 664 698 665 699 CALL init_grid_definition('boundary', grid=w_east_grid, & 666 700 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 667 701 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 668 zmin = dz, zmax = lz - dz, & 669 x0=x0, y0=y0, z0 = z0, & 670 nx = 0, ny = ny, nz = nz - 1, & 671 dx = dx, dy = dy, dz = dz) 702 x0=x0, y0=y0, z0 = z0, & 703 nx = 0, ny = ny, nz = nz - 1, z=zw) 672 704 673 705 CALL init_grid_definition('boundary', grid=w_west_grid, & 674 706 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 675 707 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 676 zmin = dz, zmax = lz - dz, & 677 x0=x0, y0=y0, z0 = z0, & 678 nx = 0, ny = ny, nz = nz - 1, & 679 dx = dx, dy = dy, dz = dz) 708 x0=x0, y0=y0, z0 = z0, & 709 nx = 0, ny = ny, nz = nz - 1, z=zw) 680 710 681 711 CALL init_grid_definition('boundary', grid=w_north_grid, & 682 712 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 683 713 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 684 zmin = dz, zmax = lz - dz, & 685 x0=x0, y0=y0, z0 = z0, & 686 nx = nx, ny = 0, nz = nz - 1, & 687 dx = dx, dy = dy, dz = dz) 714 x0=x0, y0=y0, z0 = z0, & 715 nx = nx, ny = 0, nz = nz - 1, z=zw) 688 716 689 717 CALL init_grid_definition('boundary', grid=w_south_grid, & 690 718 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 691 719 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 692 zmin = dz, zmax = lz - dz, & 693 x0=x0, y0=y0, z0 = z0, & 694 nx = nx, ny = 0, nz = nz - 1, & 695 dx = dx, dy = dy, dz = dz) 720 x0=x0, y0=y0, z0 = z0, & 721 nx = nx, ny = 0, nz = nz - 1, z=zw) 696 722 697 723 CALL init_grid_definition('boundary', grid=w_top_grid, & 698 724 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 699 725 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 700 zmin = lz, zmax = lz, & 701 x0=x0, y0=y0, z0 = z0, & 702 nx = nx, ny = ny, nz = 0, & 703 dx = dx, dy = dy, dz = dz) 726 x0=x0, y0=y0, z0 = z0, & 727 nx = nx, ny = ny, nz = 1, z=(/zw_top/)) 704 728 705 729 CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate, & 706 730 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 707 731 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 708 zmin = 0.0_dp, zmax = 0.0_dp, & 709 x0=x0, y0=y0, z0 = z0, & 710 nx = 0, ny = ny, nz = nlev - 2, & 711 dx = dx, dy = dy, dz = dz) 732 x0=x0, y0=y0, z0 = z0, & 733 nx = 0, ny = ny, nz = nlev - 2) 712 734 713 735 CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate, & 714 736 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 715 737 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 716 zmin = 0.0_dp, zmax = 0.0_dp, & 717 x0=x0, y0=y0, z0 = z0, & 718 nx = 0, ny = ny, nz = nlev - 2, & 719 dx = dx, dy = dy, dz = dz) 738 x0=x0, y0=y0, z0 = z0, & 739 nx = 0, ny = ny, nz = nlev - 2) 720 740 721 741 CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate, & 722 742 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 723 743 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 724 zmin = 0.0_dp, zmax = 0.0_dp, & 725 x0=x0, y0=y0, z0 = z0, & 726 nx = nx, ny = 0, nz = nlev - 2, & 727 dx = dx, dy = dy, dz = dz) 744 x0=x0, y0=y0, z0 = z0, & 745 nx = nx, ny = 0, nz = nlev - 2) 728 746 729 747 CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate, & 730 748 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 731 749 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 732 zmin = 0.0_dp, zmax = 0.0_dp, & 733 x0=x0, y0=y0, z0 = z0, & 734 nx = nx, ny = 0, nz = nlev - 2, & 735 dx = dx, dy = dy, dz = dz) 750 x0=x0, y0=y0, z0 = z0, & 751 nx = nx, ny = 0, nz = nlev - 2) 736 752 737 753 CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate, & 738 754 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 739 755 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 740 zmin = 0.0_dp, zmax = 0.0_dp, & 741 x0=x0, y0=y0, z0 = z0, & 742 nx = nx, ny = ny, nz = nlev - 2, & 743 dx = dx, dy = dy, dz = dz) 756 x0=x0, y0=y0, z0 = z0, & 757 nx = nx, ny = ny, nz = nlev - 2) 744 758 745 759 CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate, & 746 760 xmin = lx, xmax = lx, & 747 761 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 748 zmin = 0.0_dp, zmax = 0.0_dp, & 749 x0=x0, y0=y0, z0 = z0, & 750 nx = 0, ny = ny, nz = nlev - 2, & 751 dx = dx, dy = dy, dz = dz) 762 x0=x0, y0=y0, z0 = z0, & 763 nx = 0, ny = ny, nz = nlev - 2) 752 764 753 765 CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate, & 754 766 xmin = 0.0_dp, xmax = 0.0_dp, & 755 767 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 756 zmin = 0.0_dp, zmax = 0.0_dp, & 757 x0=x0, y0=y0, z0 = z0, & 758 nx = 0, ny = ny, nz = nlev - 2, & 759 dx = dx, dy = dy, dz = dz) 768 x0=x0, y0=y0, z0 = z0, & 769 nx = 0, ny = ny, nz = nlev - 2) 760 770 761 771 CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate, & 762 772 xmin = dx, xmax = lx - dx, & 763 773 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 764 zmin = 0.0_dp, zmax = 0.0_dp, & 765 x0=x0, y0=y0, z0 = z0, & 766 nx = nx-1, ny = 0, nz = nlev - 2, & 767 dx = dx, dy = dy, dz = dz) 774 x0=x0, y0=y0, z0 = z0, & 775 nx = nx-1, ny = 0, nz = nlev - 2) 768 776 769 777 CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate, & 770 778 xmin = dx, xmax = lx - dx, & 771 779 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 772 zmin = 0.0_dp, zmax = 0.0_dp, & 773 x0=x0, y0=y0, z0 = z0, & 774 nx = nx-1, ny = 0, nz = nlev - 2, & 775 dx = dx, dy = dy, dz = dz) 780 x0=x0, y0=y0, z0 = z0, & 781 nx = nx-1, ny = 0, nz = nlev - 2) 776 782 777 783 CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate, & 778 784 xmin = dx, xmax = lx - dx, & 779 785 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 780 zmin = 0.0_dp, zmax = 0.0_dp, & 781 x0=x0, y0=y0, z0 = z0, & 782 nx = nx-1, ny = ny, nz = nlev - 2, & 783 dx = dx, dy = dy, dz = dz) 786 x0=x0, y0=y0, z0 = z0, & 787 nx = nx-1, ny = ny, nz = nlev - 2) 784 788 785 789 CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate, & 786 790 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 787 791 ymin = dy, ymax = ly - dy, & 788 zmin = 0.0_dp, zmax = 0.0_dp, & 789 x0=x0, y0=y0, z0 = z0, & 790 nx = 0, ny = ny-1, nz = nlev - 2, & 791 dx = dx, dy = dy, dz = dz) 792 x0=x0, y0=y0, z0 = z0, & 793 nx = 0, ny = ny-1, nz = nlev - 2) 792 794 793 795 CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate, & 794 796 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 795 797 ymin = dy, ymax = ly - dy, & 796 zmin = 0.0_dp, zmax = 0.0_dp, & 797 x0=x0, y0=y0, z0 = z0, & 798 nx = 0, ny = ny-1, nz = nlev - 2, & 799 dx = dx, dy = dy, dz = dz) 798 x0=x0, y0=y0, z0 = z0, & 799 nx = 0, ny = ny-1, nz = nlev - 2) 800 800 801 801 CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate, & 802 802 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 803 803 ymin = ly, ymax = ly, & 804 zmin = 0.0_dp, zmax = 0.0_dp, & 805 x0=x0, y0=y0, z0 = z0, & 806 nx = nx, ny = 0, nz = nlev - 2, & 807 dx = dx, dy = dy, dz = dz) 804 x0=x0, y0=y0, z0 = z0, & 805 nx = nx, ny = 0, nz = nlev - 2) 808 806 809 807 CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate, & 810 808 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 811 809 ymin = 0.0_dp, ymax = 0.0_dp, & 812 zmin = 0.0_dp, zmax = 0.0_dp, & 813 x0=x0, y0=y0, z0 = z0, & 814 nx = nx, ny = 0, nz = nlev - 2, & 815 dx = dx, dy = dy, dz = dz) 810 x0=x0, y0=y0, z0 = z0, & 811 nx = nx, ny = 0, nz = nlev - 2) 816 812 817 813 CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate, & 818 814 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 819 815 ymin = dy, ymax = ly - dy, & 820 zmin = 0.0_dp, zmax = 0.0_dp, & 821 x0=x0, y0=y0, z0 = z0, & 822 nx = nx, ny = ny-1, nz = nlev - 2, & 823 dx = dx, dy = dy, dz = dz) 816 x0=x0, y0=y0, z0 = z0, & 817 nx = nx, ny = ny-1, nz = nlev - 2) 824 818 825 819 CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate, & 826 820 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 827 821 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 828 zmin = 0.0_dp, zmax = 0.0_dp, & 829 x0=x0, y0=y0, z0 = z0, & 830 nx = 0, ny = ny, nz = nlev - 2, & 831 dx = dx, dy = dy, dz = dz) 822 x0=x0, y0=y0, z0 = z0, & 823 nx = 0, ny = ny, nz = nlev - 2) 832 824 833 825 CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & 834 826 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 835 827 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 836 zmin = 0.0_dp, zmax = 0.0_dp, & 837 x0=x0, y0=y0, z0 = z0, & 838 nx = 0, ny = ny, nz = nlev - 2, & 839 dx = dx, dy = dy, dz = dz) 828 x0=x0, y0=y0, z0 = z0, & 829 nx = 0, ny = ny, nz = nlev - 2) 840 830 841 831 CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & 842 832 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 843 833 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 844 zmin = 0.0_dp, zmax = 0.0_dp, & 845 x0=x0, y0=y0, z0 = z0, & 846 nx = nx, ny = 0, nz = nlev - 2, & 847 dx = dx, dy = dy, dz = dz) 834 x0=x0, y0=y0, z0 = z0, & 835 nx = nx, ny = 0, nz = nlev - 2) 848 836 849 837 CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & 850 838 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 851 839 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 852 zmin = 0.0_dp, zmax = 0.0_dp, & 853 x0=x0, y0=y0, z0 = z0, & 854 nx = nx, ny = 0, nz = nlev - 2, & 855 dx = dx, dy = dy, dz = dz) 840 x0=x0, y0=y0, z0 = z0, & 841 nx = nx, ny = 0, nz = nlev - 2) 856 842 857 843 CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & 858 844 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 859 845 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 860 zmin = 0.0_dp, zmax = 0.0_dp, & 861 x0=x0, y0=y0, z0 = z0, & 862 nx = nx, ny = ny, nz = nlev - 2, & 863 dx = dx, dy = dy, dz = dz) 846 x0=x0, y0=y0, z0 = z0, & 847 nx = nx, ny = ny, nz = nlev - 2) 864 848 END IF 865 849 … … 869 853 !------------------------------------------------------------------------------ 870 854 871 IF (TRIM(mode) == 'profile') THEN 855 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & 856 TRIM(cfg % bc_mode) == 'ideal' ) 857 858 IF (profile_grids_required) THEN 872 859 CALL init_grid_definition('boundary', grid=scalar_profile_grid, & 873 860 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 874 861 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 875 zmin = 0.5_dp * dz, zmax = lz - 0.5_dp * dz, &876 862 x0=x0, y0=y0, z0 = z0, & 877 nx = 0, ny = 0, nz = nz, & 878 dx = dx, dy = dy, dz = dz) 863 nx = 0, ny = 0, nz = nz, z=z) 879 864 880 865 CALL init_grid_definition('boundary', grid=w_profile_grid, & 881 866 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 882 867 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 883 zmin = dz, zmax = lz - dz, &884 868 x0=x0, y0=y0, z0 = z0, & 885 nx = 0, ny = 0, nz = nz - 1, & 886 dx = dx, dy = dy, dz = dz) 869 nx = 0, ny = 0, nz = nz - 1, z=zw) 887 870 888 871 CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& 889 872 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 890 873 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 891 zmin = 0.0_dp, zmax = 0.0_dp, &892 874 x0=x0, y0=y0, z0 = z0, & 893 nx = 0, ny = 0, nz = nlev - 2, & 894 dx = dx, dy = dy, dz = dz) 875 nx = 0, ny = 0, nz = nlev - 2, z=z) 895 876 896 877 CALL init_grid_definition('boundary', grid=w_profile_intermediate, & 897 878 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 898 879 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 899 zmin = 0.0_dp, zmax = 0.0_dp, &900 880 x0=x0, y0=y0, z0 = z0, & 901 nx = 0, ny = 0, nz = nlev - 2, & 902 dx = dx, dy = dy, dz = dz) 881 nx = 0, ny = 0, nz = nlev - 2, z=zw) 903 882 END IF 904 883 … … 908 887 !------------------------------------------------------------------------------ 909 888 interp_mode = 's' 910 CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, mode=mode)889 CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode) 911 890 IF (boundary_variables_required) THEN 912 891 CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode) … … 918 897 919 898 interp_mode = 'u' 920 CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, mode=mode)899 CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 921 900 IF (boundary_variables_required) THEN 922 901 CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode) … … 928 907 929 908 interp_mode = 'v' 930 CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, mode=mode)909 CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 931 910 IF (boundary_variables_required) THEN 932 911 CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode) … … 938 917 939 918 interp_mode = 'w' 940 CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, mode=mode)919 CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 941 920 IF (boundary_variables_required) THEN 942 921 CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode) … … 947 926 END IF 948 927 949 IF (TRIM(mode) == 'profile') THEN 950 CALL setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax) 928 IF (TRIM(cfg % ic_mode) == 'profile') THEN 929 CALL setup_averaging(cosmo_grid, palm_intermediate, & 930 average_imin, average_imax, average_jmin, average_jmax) 951 931 END IF 952 932 … … 955 935 956 936 957 SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, mode)937 SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode) 958 938 959 939 TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid 960 940 TYPE(grid_definition), INTENT(INOUT), TARGET :: palm_grid, palm_intermediate 961 941 CHARACTER, INTENT(IN) :: kind 962 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: mode942 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 963 943 964 944 TYPE(grid_definition), POINTER :: grid … … 966 946 REAL(dp), DIMENSION(:,:,:), POINTER :: h 967 947 968 LOGICAL :: setup_vertical = .TRUE. 969 970 IF (PRESENT(mode)) THEN 971 IF (TRIM(mode) == 'profile') setup_vertical = .FALSE. 972 ELSE 973 setup_vertical = .TRUE. 948 LOGICAL :: setup_vertical 949 950 setup_vertical = .TRUE. 951 IF (PRESENT(ic_mode)) THEN 952 IF (TRIM(ic_mode) == 'profile') setup_vertical = .FALSE. 974 953 END IF 975 954 … … 1006 985 CASE DEFAULT 1007 986 1008 message = "Interpolation mode '" // mode// "' is not supported."987 message = "Interpolation quantity '" // kind // "' is not supported." 1009 988 CALL abort('setup_interpolation', message) 1010 989 … … 1014 993 1015 994 CALL find_horizontal_neighbours(lat, lon, & 1016 cosmo_grid % dxi, cosmo_grid % dyi, grid % clat, & 1017 grid % clon, grid % ii, grid % jj) 995 grid % clat, grid % clon, grid % ii, grid % jj) 1018 996 1019 997 CALL compute_horizontal_interp_weights(lat, lon, & 1020 cosmo_grid % dxi, cosmo_grid % dyi, grid % clat, & 1021 grid % clon, grid % ii, grid % jj, grid % w_horiz) 998 grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz) 1022 999 1023 1000 !------------------------------------------------------------------------------ … … 1044 1021 1045 1022 TYPE(grid_definition), POINTER :: grid 1046 REAL :: lonmin_pos,lonmax_pos, latmin_pos, latmax_pos 1023 REAL(dp) :: lonmin_pos,lonmax_pos, latmin_pos, latmax_pos 1024 REAL(dp) :: cosmo_dxi, cosmo_dyi 1025 1026 cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0)) 1027 cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0)) 1047 1028 1048 1029 ! find horizontal index ranges for profile averaging 1049 lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_ grid %dxi1050 lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_ grid %dxi1051 latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_ grid %dyi1052 latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_ grid %dyi1030 lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi 1031 lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi 1032 latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi 1033 latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi 1053 1034 1054 1035 imin = FLOOR(lonmin_pos) … … 1076 1057 ! Description: 1077 1058 ! ------------ 1078 !> Helper function that computes horizontal domain extend in x or y direction1079 !> such that the centres of a boundary grid fall at -dx/2 or lx + dx/2.1080 !>1081 !> Input parameters:1082 !> -----------------1083 !> dxy : grid spacing in x or y direction1084 !> lxy : domain length in dxy direction1085 !>1086 !> Output parameters:1087 !> ------------------1088 !> boundary_extent : Domain minimum xymin (maximum xymax) if dxy < 0 (> 0)1089 !------------------------------------------------------------------------------!1090 REAL(dp) FUNCTION boundary_extent(dxy, lxy)1091 REAL(dp), INTENT(IN) :: dxy, lxy1092 1093 boundary_extent = 0.5_dp * lxy + SIGN(lxy + ABS(dxy), dxy)1094 1095 END FUNCTION boundary_extent1096 1097 1098 !------------------------------------------------------------------------------!1099 ! Description:1100 ! ------------1101 1059 !> Initializes grid_definition-type variables. 1102 1060 !> 1103 1061 !> Input parameters: 1104 1062 !> ----------------- 1105 !> mode : Initialization mode, distinguishes between PALM-4U and COSMO-DE grids1106 !> as well as grids covering the boundary surfaces. Valid modes are:1063 !> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids 1064 !> as well as grids covering the boundary surfaces. Valid kinds are: 1107 1065 !> - 'palm' 1108 1066 !> - 'cosmo-de' … … 1114 1072 !> PALM-4U computational domain (i.e. the outer cell faces). The coordinates 1115 1073 !> of the generated grid will be inferred from this information taking into 1116 !> account the initialization mode . For example, the coordinates of a1074 !> account the initialization mode ic_mode. For example, the coordinates of a 1117 1075 !> boundary grid initialized using mode 'eastwest-scalar' will be located in 1118 1076 !> planes one half grid point outwards of xmin and xmax. … … 1126 1084 !> grid : Grid variable to be initialized. 1127 1085 !------------------------------------------------------------------------------! 1128 SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, zmin, zmax,&1129 x0, y0, z0, nx, ny, nz, dx, dy, dz, grid,mode)1086 SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, & 1087 x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode) 1130 1088 CHARACTER(LEN=*), INTENT(IN) :: kind 1131 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: mode1089 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 1132 1090 INTEGER, INTENT(IN) :: nx, ny, nz 1133 REAL(dp), INTENT(IN) :: xmin, xmax, ymin, ymax , zmin, zmax1091 REAL(dp), INTENT(IN) :: xmin, xmax, ymin, ymax 1134 1092 REAL(dp), INTENT(IN) :: x0, y0, z0 1135 REAL(dp), OPTIONAL, INTENT(IN) :: dx, dy, dz 1093 REAL(dp), INTENT(IN), TARGET, OPTIONAL :: z(:) 1094 REAL(dp), INTENT(IN), TARGET, OPTIONAL :: zw(:) 1136 1095 TYPE(grid_definition), INTENT(INOUT) :: grid 1137 1096 … … 1142 1101 grid % lx = xmax - xmin 1143 1102 grid % ly = ymax - ymin 1144 grid % lz = zmax - zmin1145 1103 1146 1104 grid % x0 = x0 … … 1151 1109 1152 1110 CASE('boundary') 1153 IF (.NOT. PRESENT(dx)) THEN 1154 message = "dx is not present but needed for 'eastwest-scalar' "//&1155 "initializaton."1111 1112 IF (.NOT.PRESENT(z)) THEN 1113 message = "z has not been passed but is required for 'boundary' grids" 1156 1114 CALL abort('init_grid_definition', message) 1157 1115 END IF 1158 IF (.NOT. PRESENT(dy)) THEN1159 message = "dy is not present but needed for 'eastwest-scalar' "//&1160 "initializaton."1161 CALL abort('init_grid_definition', message)1162 END IF1163 IF (.NOT. PRESENT(dz)) THEN1164 message = "dz is not present but needed for 'eastwest-scalar' "//&1165 "initializaton."1166 CALL abort('init_grid_definition', message)1167 END IF1168 1169 grid % dx = dx1170 grid % dy = dy1171 grid % dz = dz1172 1173 grid % dxi = 1.0_dp / grid % dx1174 grid % dyi = 1.0_dp / grid % dy1175 grid % dzi = 1.0_dp / grid % dz1176 1116 1177 1117 ALLOCATE( grid % x(0:nx) ) … … 1181 1121 CALL linspace(ymin, ymax, grid % y) 1182 1122 1183 ALLOCATE( grid % z(0:nz) ) 1184 CALL linspace(zmin, zmax, grid % z) 1123 grid % z => z 1185 1124 1186 1125 ! Allocate neighbour indices and weights 1187 IF (TRIM( mode) .NE. 'profile') THEN1188 ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )1126 IF (TRIM(ic_mode) .NE. 'profile') THEN 1127 ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) 1189 1128 grid % kk(:,:,:,:) = -1 1190 1129 1191 ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )1130 ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) ) 1192 1131 grid % w_verti(:,:,:,:) = 0.0_dp 1193 1132 END IF 1194 1133 1195 1134 CASE('boundary intermediate') 1196 IF (.NOT. PRESENT(dx)) THEN1197 message = "dx is not present but needed for 'eastwest-scalar' "//&1198 "initializaton."1199 CALL abort('init_grid_definition', message)1200 END IF1201 IF (.NOT. PRESENT(dy)) THEN1202 message = "dy is not present but needed for 'eastwest-scalar' "//&1203 "initializaton."1204 CALL abort('init_grid_definition', message)1205 END IF1206 IF (.NOT. PRESENT(dz)) THEN1207 message = "dz is not present but needed for 'eastwest-scalar' "//&1208 "initializaton."1209 CALL abort('init_grid_definition', message)1210 END IF1211 1212 grid % dx = dx1213 grid % dy = dy1214 grid % dz = dz1215 1216 grid % dxi = 1.0_dp / grid % dx1217 grid % dyi = 1.0_dp / grid % dy1218 grid % dzi = 1.0_dp / grid % dz1219 1135 1220 1136 ALLOCATE( grid % x(0:nx) ) … … 1224 1140 CALL linspace(ymin, ymax, grid % y) 1225 1141 1226 ALLOCATE( grid % z(0:nz) )1227 CALL linspace(zmin, zmax, grid % z)1228 1229 1142 ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) ) 1143 1230 1144 CALL rotate_to_cosmo( & 1231 phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude1232 lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude1145 phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude 1146 lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude 1233 1147 phip = phi_cn, lamp = lambda_cn, & 1234 1148 phi = grid % clat, & … … 1249 1163 ! corresponding latitudes and longitudes of the rotated pole grid. 1250 1164 CASE('palm') 1165 1166 IF (.NOT.PRESENT(z)) THEN 1167 message = "z has not been passed but is required for 'palm' grids" 1168 CALL abort('init_grid_definition', message) 1169 END IF 1170 1171 IF (.NOT.PRESENT(zw)) THEN 1172 message = "zw has not been passed but is required for 'palm' grids" 1173 CALL abort('init_grid_definition', message) 1174 END IF 1175 1251 1176 grid % name(1) = 'x and lon' 1252 1177 grid % name(2) = 'y and lat' 1253 1178 grid % name(3) = 'z' 1254 1179 1255 grid % dx = grid % lx / (nx + 1) 1256 grid % dy = grid % ly / (ny + 1) 1257 grid % dz = grid % lz / (nz + 1) 1258 1259 grid % dxi = 1.0_dp / grid % dx 1260 grid % dyi = 1.0_dp / grid % dy 1261 grid % dzi = 1.0_dp / grid % dz 1262 1263 ALLOCATE( grid % x(0:nx), grid % y(0:ny), grid % z(0:nz) ) 1264 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny), grid % zw(1:nz) ) 1265 CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x) 1266 CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y) 1267 CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z) 1268 CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu) 1269 CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv) 1270 CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw) 1180 !TODO: Remove use of global dx, dy, dz variables. Consider 1181 !TODO: associating global x,y, and z arrays. 1182 ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) 1183 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) 1184 CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x) 1185 CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y) 1186 grid % z => z 1187 CALL linspace(xmin + dx, xmax - dx, grid % xu) 1188 CALL linspace(ymin + dy, ymax - dy, grid % yv) 1189 grid % zw => zw 1271 1190 1272 1191 grid % depths => depths 1273 1192 1274 1193 ! Allocate neighbour indices and weights 1275 IF (TRIM( mode) .NE. 'profile') THEN1276 ALLOCATE( grid % kk(0:nx, 0:ny, 0:nz, 2) )1194 IF (TRIM(ic_mode) .NE. 'profile') THEN 1195 ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) 1277 1196 grid % kk(:,:,:,:) = -1 1278 1197 1279 ALLOCATE( grid % w_verti(0:nx, 0:ny, 0:nz, 2) )1198 ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) ) 1280 1199 grid % w_verti(:,:,:,:) = 0.0_dp 1281 1200 END IF 1282 1201 1283 1202 CASE('palm intermediate') 1203 1284 1204 grid % name(1) = 'x and lon' 1285 1205 grid % name(2) = 'y and lat' 1286 grid % name(3) = 'z' 1287 1288 grid % dx = grid % lx / (nx + 1) 1289 grid % dy = grid % ly / (ny + 1) 1290 grid % dz = grid % lz / (nz + 1) 1291 1292 grid % dxi = 1.0_dp / grid % dx 1293 grid % dyi = 1.0_dp / grid % dy 1294 grid % dzi = 1.0_dp / grid % dz 1295 1296 ALLOCATE( grid % x(0:nx), grid % y(0:ny), grid % z(0:nz) ) 1297 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny), grid % zw(1:nz) ) 1298 CALL linspace(xmin + 0.5_dp*grid % dx, xmax - 0.5_dp*grid % dx, grid % x) 1299 CALL linspace(ymin + 0.5_dp*grid % dy, ymax - 0.5_dp*grid % dy, grid % y) 1300 CALL linspace(zmin + 0.5_dp*grid % dz, zmax - 0.5_dp*grid % dz, grid % z) 1301 CALL linspace(xmin + grid % dx, xmax - grid % dx, grid % xu) 1302 CALL linspace(ymin + grid % dy, ymax - grid % dy, grid % yv) 1303 CALL linspace(zmin + grid % dz, zmax - grid % dz, grid % zw) 1206 grid % name(3) = 'interpolated hhl or hfl' 1207 1208 !TODO: Remove use of global dx, dy, dz variables. Consider 1209 !TODO: associating global x,y, and z arrays. 1210 ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) 1211 ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) 1212 CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid % x) 1213 CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid % y) 1214 CALL linspace(xmin + dx, xmax - dx, grid % xu) 1215 CALL linspace(ymin + dy, ymax - dy, grid % yv) 1304 1216 1305 1217 grid % depths => depths … … 1355 1267 grid % name(3) = 'height' 1356 1268 1357 grid % dx = grid % lx / nx ! = 0.025 deg, stored in radians 1358 grid % dy = grid % ly / ny ! = 0.025 deg, stored in radians 1359 grid % dz = 0.0_dp ! not defined yet 1360 1361 grid % dxi = 1.0_dp / grid % dx ! [rad^-1] 1362 grid % dyi = 1.0_dp / grid % dy ! [rad^-1] 1363 grid % dzi = 0.0_dp ! not defined yet 1364 1365 ALLOCATE( grid % lon(0:nx), grid % lat(0:ny), grid % z(0:nz) ) 1366 ALLOCATE( grid % lonu(0:nx), grid % latv(0:ny), grid % zw(0:nz) ) 1269 ALLOCATE( grid % lon(0:nx), grid % lat(0:ny) ) 1270 ALLOCATE( grid % lonu(0:nx), grid % latv(0:ny) ) 1367 1271 1368 1272 CALL linspace(xmin, xmax, grid % lon) 1369 1273 CALL linspace(ymin, ymax, grid % lat) 1370 grid % lonu(:) = grid % lon + 0.5_dp * grid % dx1371 grid % latv(:) = grid % lat + 0.5_dp * grid % dy1274 grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx) 1275 grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny) 1372 1276 1373 1277 ! Point to heights of half levels (hhl) and compute heights of full … … 1384 1288 1385 1289 END SUBROUTINE init_grid_definition 1290 1291 1292 !------------------------------------------------------------------------------! 1293 ! Description: 1294 ! ------------ 1295 !> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see 1296 !> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139 1297 !> 1298 !> This routine computes the levels of scalar points. The levels of the velocity 1299 !> points are then obtained as the midpoints inbetween using the INIFOR routine 1300 !> 'modpoints'. 1301 !------------------------------------------------------------------------------! 1302 SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, & 1303 dz_stretch_level_start, dz_stretch_level_end, & 1304 dz_stretch_factor_array) 1305 1306 REAL(dp), DIMENSION(:), INTENT(INOUT) :: z, dz, dz_stretch_factor_array 1307 REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_level_start, dz_stretch_level_end 1308 REAL(dp), INTENT(IN) :: dz_max, dz_stretch_factor, dz_stretch_level 1309 1310 INTEGER :: number_stretch_level_start !< number of user-specified start levels for stretching 1311 INTEGER :: number_stretch_level_end !< number of user-specified end levels for stretching 1312 1313 REAL(dp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end 1314 REAL(dp) :: dz_level_end, dz_stretched 1315 1316 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 1317 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 1318 INTEGER :: dz_stretch_level_index = 0 1319 INTEGER :: k, n, number_dz 1320 ! 1321 !-- Compute height of u-levels from constant grid length and dz stretch factors 1322 IF ( dz(1) == -1.0_dp ) THEN 1323 message = 'missing dz' 1324 CALL abort( 'stretched_z', message) 1325 ELSEIF ( dz(1) <= 0.0_dp ) THEN 1326 WRITE( message, * ) 'dz=', dz(1),' <= 0.0' 1327 CALL abort( 'stretched_z', message) 1328 ENDIF 1329 1330 ! 1331 !-- Initialize dz_stretch_level_start with the value of dz_stretch_level 1332 !-- if it was set by the user 1333 IF ( dz_stretch_level /= -9999999.9_dp ) THEN 1334 dz_stretch_level_start(1) = dz_stretch_level 1335 ENDIF 1336 1337 ! 1338 !-- Determine number of dz values and stretching levels specified by the 1339 !-- user to allow right controlling of the stretching mechanism and to 1340 !-- perform error checks. The additional requirement that dz /= dz_max 1341 !-- for counting number of user-specified dz values is necessary. Otherwise 1342 !-- restarts would abort if the old stretching mechanism with dz_stretch_level 1343 !-- is used (Attention: The user is not allowed to specify a dz value equal 1344 !-- to the default of dz_max = 999.0). 1345 number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max ) 1346 number_stretch_level_start = COUNT( dz_stretch_level_start /= & 1347 -9999999.9_dp ) 1348 number_stretch_level_end = COUNT( dz_stretch_level_end /= & 1349 9999999.9_dp ) 1350 1351 ! 1352 !-- The number of specified end levels +1 has to be the same than the number 1353 !-- of specified dz values 1354 IF ( number_dz /= number_stretch_level_end + 1 ) THEN 1355 WRITE( message, * ) 'The number of values for dz = ', & 1356 number_dz, 'has to be the same than ', & 1357 'the number of values for ', & 1358 'dz_stretch_level_end + 1 = ', & 1359 number_stretch_level_end+1 1360 CALL abort( 'stretched_z', message) 1361 ENDIF 1362 1363 ! 1364 !-- The number of specified start levels has to be the same or one less than 1365 !-- the number of specified dz values 1366 IF ( number_dz /= number_stretch_level_start + 1 .AND. & 1367 number_dz /= number_stretch_level_start ) THEN 1368 WRITE( message, * ) 'The number of values for dz = ', & 1369 number_dz, 'has to be the same or one ', & 1370 'more than& the number of values for ', & 1371 'dz_stretch_level_start = ', & 1372 number_stretch_level_start 1373 CALL abort( 'stretched_z', message) 1374 ENDIF 1375 1376 !-- The number of specified start levels has to be the same or one more than 1377 !-- the number of specified end levels 1378 IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. & 1379 number_stretch_level_start /= number_stretch_level_end ) THEN 1380 WRITE( message, * ) 'The number of values for ', & 1381 'dz_stretch_level_start = ', & 1382 dz_stretch_level_start, 'has to be the ',& 1383 'same or one more than& the number of ', & 1384 'values for dz_stretch_level_end = ', & 1385 number_stretch_level_end 1386 CALL abort( 'stretched_z', message) 1387 ENDIF 1388 1389 ! 1390 !-- Initialize dz for the free atmosphere with the value of dz_max 1391 IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND. & 1392 number_stretch_level_start /= 0 ) THEN 1393 dz(number_stretch_level_start+1) = dz_max 1394 ENDIF 1395 1396 ! 1397 !-- Initialize the stretching factor if (infinitely) stretching in the free 1398 !-- atmosphere is desired (dz_stretch_level_end was not specified for the 1399 !-- free atmosphere) 1400 IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 1401 dz_stretch_factor_array(number_stretch_level_start) = & 1402 dz_stretch_factor 1403 ENDIF 1404 1405 !-- Allocation of arrays for stretching 1406 ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) ) 1407 1408 ! 1409 !-- The stretching region has to be large enough to allow for a smooth 1410 !-- transition between two different grid spacings 1411 DO n = 1, number_stretch_level_start 1412 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & 1413 4 * MAX( dz(n),dz(n+1) ) 1414 ENDDO 1415 1416 IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > & 1417 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN 1418 !IF ( ANY( min_dz_stretch_level_end > & 1419 ! dz_stretch_level_end ) ) THEN 1420 message = 'Each dz_stretch_level_end has to be larger ' // & 1421 'than its corresponding value for ' // & 1422 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& 1423 'to allow for smooth grid stretching' 1424 CALL abort('stretched_z', message) 1425 ENDIF 1426 1427 ! 1428 !-- Stretching must not be applied within the prandtl_layer 1429 !-- (first two grid points). For the default case dz_stretch_level_start 1430 !-- is negative. Therefore the absolut value is checked here. 1431 IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN 1432 WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',& 1433 'larger than ', dz(1) * 1.5 1434 CALL abort( 'stretched_z', message) 1435 ENDIF 1436 1437 ! 1438 !-- The stretching has to start and end on a grid level. Therefore 1439 !-- user-specified values have to ''interpolate'' to the next lowest level 1440 IF ( number_stretch_level_start /= 0 ) THEN 1441 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & 1442 dz(1)/2.0) / dz(1) ) & 1443 * dz(1) + dz(1)/2.0 1444 ENDIF 1445 1446 IF ( number_stretch_level_start > 1 ) THEN 1447 DO n = 2, number_stretch_level_start 1448 dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & 1449 dz(n) ) * dz(n) 1450 ENDDO 1451 ENDIF 1452 1453 IF ( number_stretch_level_end /= 0 ) THEN 1454 DO n = 1, number_stretch_level_end 1455 dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & 1456 dz(n+1) ) * dz(n+1) 1457 ENDDO 1458 ENDIF 1459 1460 ! 1461 !-- Determine stretching factor if necessary 1462 IF ( number_stretch_level_end >= 1 ) THEN 1463 CALL calculate_stretching_factor( number_stretch_level_end, dz, & 1464 dz_stretch_factor, & 1465 dz_stretch_factor_array, & 1466 dz_stretch_level_end, & 1467 dz_stretch_level_start ) 1468 ENDIF 1469 1470 z(1) = dz(1) * 0.5_dp 1471 ! 1472 dz_stretch_level_index = n 1473 dz_stretched = dz(1) 1474 DO k = 2, n 1475 1476 IF ( dz_stretch_level <= z(k-1) .AND. dz_stretched < dz_max ) THEN 1477 1478 dz_stretched = dz_stretched * dz_stretch_factor 1479 dz_stretched = MIN( dz_stretched, dz_max ) 1480 1481 IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1 1482 1483 ENDIF 1484 1485 z(k) = z(k-1) + dz_stretched 1486 1487 ENDDO 1488 !-- Determine u and v height levels considering the possibility of grid 1489 !-- stretching in several heights. 1490 n = 1 1491 dz_stretch_level_start_index(:) = UBOUND(z, 1) 1492 dz_stretch_level_end_index(:) = UBOUND(z, 1) 1493 dz_stretched = dz(1) 1494 1495 !-- The default value of dz_stretch_level_start is negative, thus the first 1496 !-- condition is always true. Hence, the second condition is necessary. 1497 DO k = 2, UBOUND(z, 1) 1498 IF ( dz_stretch_level_start(n) <= z(k-1) .AND. & 1499 dz_stretch_level_start(n) /= -9999999.9_dp ) THEN 1500 dz_stretched = dz_stretched * dz_stretch_factor_array(n) 1501 1502 IF ( dz(n) > dz(n+1) ) THEN 1503 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz 1504 ELSE 1505 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz 1506 ENDIF 1507 1508 IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) ) & 1509 dz_stretch_level_start_index(n) = k-1 1510 1511 ENDIF 1512 1513 z(k) = z(k-1) + dz_stretched 1514 1515 ! 1516 !-- Make sure that the stretching ends exactly at dz_stretch_level_end 1517 dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) 1518 1519 IF ( dz_level_end < dz(n+1)/3.0 ) THEN 1520 z(k) = dz_stretch_level_end(n) 1521 dz_stretched = dz(n+1) 1522 dz_stretch_level_end_index(n) = k 1523 n = n + 1 1524 ENDIF 1525 ENDDO 1526 1527 DEALLOCATE( min_dz_stretch_level_end ) 1528 1529 END SUBROUTINE stretched_z 1530 1531 1532 ! Description: [PALM subroutine] 1533 ! -----------------------------------------------------------------------------! 1534 !> Calculation of the stretching factor through an iterative method. Ideas were 1535 !> taken from the paper "Regional stretched grid generation and its application 1536 !> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the 1537 !> system of equations has two variables (r,l) but four requirements 1538 !> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which 1539 !> results into an overdetermined system. 1540 !------------------------------------------------------------------------------! 1541 SUBROUTINE calculate_stretching_factor( number_end, dz, dz_stretch_factor, & 1542 dz_stretch_factor_array, & 1543 dz_stretch_level_end, & 1544 dz_stretch_level_start ) 1545 1546 REAL(dp), DIMENSION(:), INTENT(IN) :: dz 1547 REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_factor_array 1548 REAL(dp), DIMENSION(:), INTENT(IN) :: dz_stretch_level_end, dz_stretch_level_start 1549 REAL(dp) :: dz_stretch_factor 1550 1551 INTEGER :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached 1552 INTEGER :: l_rounded !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2 1553 INTEGER :: n !< loop variable for stretching 1554 1555 INTEGER, INTENT(IN) :: number_end !< number of user-specified end levels for stretching 1556 1557 REAL(dp) :: delta_l !< absolute difference between l and l_rounded 1558 REAL(dp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2 1559 REAL(dp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) 1560 REAL(dp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration 1561 REAL(dp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region) 1562 REAL(dp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly 1563 REAL(dp) :: numerator !< numerator of the quotient 1564 REAL(dp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly 1565 REAL(dp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly 1566 1567 REAL(dp) :: dz_stretch_factor_array_2(9) = 1.08_dp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 1568 1569 REAL(dp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors 1570 REAL(dp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor 1571 REAL(dp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor 1572 1573 1574 l = 0 1575 DO n = 1, number_end 1576 1577 iterations = 1 1578 stretch_factor_1 = 1.0 1579 stretch_factor_2 = 1.0 1580 delta_total_old = 1.0 1581 1582 IF ( dz(n) > dz(n+1) ) THEN 1583 DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 1584 1585 stretch_factor_1 = 1.0 - iterations * stretch_factor_interval 1586 distance = ABS( dz_stretch_level_end(n) - & 1587 dz_stretch_level_start(n) ) 1588 numerator = distance*stretch_factor_1/dz(n) + & 1589 stretch_factor_1 - distance/dz(n) 1590 1591 IF ( numerator > 0.0 ) THEN 1592 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1593 l_rounded = NINT( l ) 1594 delta_l = ABS( l_rounded - l ) / l 1595 ENDIF 1596 1597 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1598 1599 delta_stretch_factor = ABS( stretch_factor_1 - & 1600 stretch_factor_2 ) / & 1601 stretch_factor_2 1602 1603 delta_total_new = delta_l + delta_stretch_factor 1604 1605 ! 1606 !-- stretch_factor_1 is taken to guarantee that the stretching 1607 !-- procedure ends as close as possible to dz_stretch_level_end. 1608 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1609 !-- equal to dz(n+1) after l_rounded grid levels. 1610 IF (delta_total_new < delta_total_old) THEN 1611 dz_stretch_factor_array(n) = stretch_factor_1 1612 dz_stretch_factor_array_2(n) = stretch_factor_2 1613 delta_total_old = delta_total_new 1614 ENDIF 1615 1616 iterations = iterations + 1 1617 1618 ENDDO 1619 1620 ELSEIF ( dz(n) < dz(n+1) ) THEN 1621 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) 1622 1623 stretch_factor_1 = 1.0 + iterations * stretch_factor_interval 1624 distance = ABS( dz_stretch_level_end(n) - & 1625 dz_stretch_level_start(n) ) 1626 numerator = distance*stretch_factor_1/dz(n) + & 1627 stretch_factor_1 - distance/dz(n) 1628 1629 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1630 l_rounded = NINT( l ) 1631 delta_l = ABS( l_rounded - l ) / l 1632 1633 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1634 1635 delta_stretch_factor = ABS( stretch_factor_1 - & 1636 stretch_factor_2 ) / & 1637 stretch_factor_2 1638 1639 delta_total_new = delta_l + delta_stretch_factor 1640 1641 ! 1642 !-- stretch_factor_1 is taken to guarantee that the stretching 1643 !-- procedure ends as close as possible to dz_stretch_level_end. 1644 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1645 !-- equal to dz(n+1) after l_rounded grid levels. 1646 IF (delta_total_new < delta_total_old) THEN 1647 dz_stretch_factor_array(n) = stretch_factor_1 1648 dz_stretch_factor_array_2(n) = stretch_factor_2 1649 delta_total_old = delta_total_new 1650 ENDIF 1651 1652 iterations = iterations + 1 1653 ENDDO 1654 1655 ELSE 1656 message = 'Two adjacent values of dz must be different' 1657 CALL abort( 'calculate_stretching_factor', message) 1658 ENDIF 1659 1660 ! 1661 !-- Check if also the second stretching factor fits into the allowed 1662 !-- interval. If not, print a warning for the user. 1663 IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR. & 1664 dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN 1665 WRITE( message, * ) 'stretch_factor_2 = ', & 1666 dz_stretch_factor_array_2(n), ' which is',& 1667 ' responsible for exactly reaching& dz =',& 1668 dz(n+1), 'after a specific amount of', & 1669 ' grid levels& exceeds the upper', & 1670 ' limit =', stretch_factor_upper_limit, & 1671 ' &or lower limit = ', & 1672 stretch_factor_lower_limit 1673 CALL abort( 'calculate_stretching_factor', message ) 1674 1675 ENDIF 1676 ENDDO 1677 1678 END SUBROUTINE calculate_stretching_factor 1679 1680 SUBROUTINE midpoints(z, zw) 1681 1682 REAL(dp), INTENT(IN) :: z(0:) 1683 REAL(dp), INTENT(OUT) :: zw(1:) 1684 1685 INTEGER :: k 1686 1687 DO k = 1, UBOUND(zw, 1) 1688 zw(k) = 0.5_dp * (z(k-1) + z(k)) 1689 END DO 1690 1691 END SUBROUTINE midpoints 1386 1692 1387 1693 … … 1409 1715 ) 1410 1716 1411 !potential temperature, surface pressure 1717 !potential temperature, surface pressure, including nudging and subsidence 1412 1718 io_group_list(3) = init_io_group( & 1413 1719 in_files = flow_files, & 1414 out_vars = [output_var_table(3:8), output_var_table(42:42)], & 1720 out_vars = [output_var_table(3:8), output_var_table(42:42), & 1721 output_var_table(49:51)], & 1415 1722 in_var_list = (/input_var_table(3), input_var_table(17)/), & 1416 1723 kind = 'temperature' & 1417 1724 ) 1418 1725 1419 !specific humidity 1726 !specific humidity including nudging and subsidence 1420 1727 io_group_list(4) = init_io_group( & 1421 1728 in_files = flow_files, & 1422 out_vars = output_var_table(9:14),&1729 out_vars = [output_var_table(9:14), output_var_table(52:54)], & 1423 1730 in_var_list = input_var_table(4:4), & 1424 1731 kind = 'scalar' & … … 1428 1735 io_group_list(5) = init_io_group( & 1429 1736 in_files = flow_files, & 1430 out_vars = [output_var_table(15:26), output_var_table(43:4 4)], &1737 out_vars = [output_var_table(15:26), output_var_table(43:46)], & 1431 1738 !out_vars = output_var_table(15:20), & 1432 1739 in_var_list = input_var_table(5:6), & … … 1444 1751 !io_group_list(6) % to_be_processed = .FALSE. 1445 1752 1446 !w velocity 1753 !w velocity and subsidence and w nudging 1447 1754 io_group_list(7) = init_io_group( & 1448 1755 in_files = flow_files, & 1449 out_vars = output_var_table(27:32),&1756 out_vars = [output_var_table(27:32), output_var_table(47:48)], & 1450 1757 in_var_list = input_var_table(7:7), & 1451 1758 kind = 'scalar' & … … 1459 1766 kind = 'accumulated' & 1460 1767 ) 1768 io_group_list(8) % to_be_processed = .FALSE. 1461 1769 1462 1770 !snow … … 1467 1775 kind = 'accumulated' & 1468 1776 ) 1777 io_group_list(9) % to_be_processed = .FALSE. 1469 1778 1470 1779 !graupel … … 1475 1784 kind = 'accumulated' & 1476 1785 ) 1786 io_group_list(10) % to_be_processed = .FALSE. 1477 1787 1478 1788 !evapotranspiration … … 1483 1793 kind = 'accumulated' & 1484 1794 ) 1795 io_group_list(11) % to_be_processed = .FALSE. 1485 1796 1486 1797 !2m air temperature … … 1517 1828 kind = 'running average' & 1518 1829 ) 1830 io_group_list(15) % to_be_processed = .FALSE. 1519 1831 1520 1832 !lw radiation balance … … 1525 1837 kind = 'running average' & 1526 1838 ) 1839 io_group_list(16) % to_be_processed = .FALSE. 1527 1840 1528 1841 END SUBROUTINE setup_io_groups … … 1562 1875 CALL report('fini_grids', 'Deallocating grids') 1563 1876 1877 DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column) 1878 1564 1879 DEALLOCATE(palm_grid%x, palm_grid%y, palm_grid%z, & 1565 1880 palm_grid%xu, palm_grid%yv, palm_grid%zw, & … … 1572 1887 palm_intermediate%clonu, palm_intermediate%clatu) 1573 1888 1574 DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat, cosmo_grid%z,&1575 cosmo_grid%lonu, cosmo_grid%latv, cosmo_grid%zw,&1889 DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat, & 1890 cosmo_grid%lonu, cosmo_grid%latv, & 1576 1891 cosmo_grid%hfl) 1577 1892 … … 1584 1899 !> Initializes the the variable list. 1585 1900 !------------------------------------------------------------------------------! 1586 SUBROUTINE setup_variable_tables( mode)1587 CHARACTER(LEN=*), INTENT(IN) :: mode1901 SUBROUTINE setup_variable_tables(ic_mode) 1902 CHARACTER(LEN=*), INTENT(IN) :: ic_mode 1588 1903 TYPE(nc_var), POINTER :: var 1589 1904 1590 IF (TRIM( start_date) == '') THEN1905 IF (TRIM(cfg % start_date) == '') THEN 1591 1906 message = 'Simulation start date has not been set.' 1592 1907 CALL abort('setup_variable_tables', message) 1593 1908 END IF 1594 1909 1595 nc_source_text = 'COSMO-DE analysis from ' // TRIM( start_date)1910 nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date) 1596 1911 1597 1912 n_invar = 17 1598 n_outvar = 441913 n_outvar = 55 1599 1914 ALLOCATE( input_var_table(n_invar) ) 1600 1915 ALLOCATE( output_var_table(n_outvar) ) … … 1693 2008 !- Section 2: NetCDF output variables 1694 2009 !------------------------------------------------------------------------------ 2010 ! 2011 !------------------------------------------------------------------------------ 2012 ! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions 2013 !------------------------------------------------------------------------------ 1695 2014 output_var_table(1) = init_nc_var( & 1696 2015 name = 'init_soil_t', & … … 1718 2037 1719 2038 output_var_table(3) = init_nc_var( & 1720 name = 'init_ pt',&2039 name = 'init_atmosphere_pt', & 1721 2040 std_name = "", & 1722 2041 long_name = "initial potential temperature", & … … 1727 2046 grid = palm_grid, & 1728 2047 intermediate_grid = palm_intermediate, & 1729 is_profile = (TRIM( mode) == 'profile')&1730 ) 1731 IF (TRIM( mode) == 'profile') THEN2048 is_profile = (TRIM(ic_mode) == 'profile') & 2049 ) 2050 IF (TRIM(ic_mode) == 'profile') THEN 1732 2051 output_var_table(3) % grid => scalar_profile_grid 1733 2052 output_var_table(3) % intermediate_grid => scalar_profile_intermediate … … 1795 2114 1796 2115 output_var_table(9) = init_nc_var( & 1797 name = 'init_ qv',&2116 name = 'init_atmosphere_qv', & 1798 2117 std_name = "", & 1799 2118 long_name = "initial specific humidity", & … … 1804 2123 grid = palm_grid, & 1805 2124 intermediate_grid = palm_intermediate, & 1806 is_profile = (TRIM( mode) == 'profile')&1807 ) 1808 IF (TRIM( mode) == 'profile') THEN2125 is_profile = (TRIM(ic_mode) == 'profile') & 2126 ) 2127 IF (TRIM(ic_mode) == 'profile') THEN 1809 2128 output_var_table(9) % grid => scalar_profile_grid 1810 2129 output_var_table(9) % intermediate_grid => scalar_profile_intermediate … … 1872 2191 1873 2192 output_var_table(15) = init_nc_var( & 1874 name = 'init_ u',&2193 name = 'init_atmosphere_u', & 1875 2194 std_name = "", & 1876 2195 long_name = "initial wind component in x direction", & … … 1881 2200 grid = u_initial_grid, & 1882 2201 intermediate_grid = u_initial_intermediate, & 1883 is_profile = (TRIM( mode) == 'profile')&1884 ) 1885 IF (TRIM( mode) == 'profile') THEN2202 is_profile = (TRIM(ic_mode) == 'profile') & 2203 ) 2204 IF (TRIM(ic_mode) == 'profile') THEN 1886 2205 output_var_table(15) % grid => scalar_profile_grid 1887 2206 output_var_table(15) % intermediate_grid => scalar_profile_intermediate … … 1949 2268 1950 2269 output_var_table(21) = init_nc_var( & 1951 name = 'init_ v',&2270 name = 'init_atmosphere_v', & 1952 2271 std_name = "", & 1953 2272 long_name = "initial wind component in y direction", & … … 1958 2277 grid = v_initial_grid, & 1959 2278 intermediate_grid = v_initial_intermediate, & 1960 is_profile = (TRIM( mode) == 'profile')&1961 ) 1962 IF (TRIM( mode) == 'profile') THEN2279 is_profile = (TRIM(ic_mode) == 'profile') & 2280 ) 2281 IF (TRIM(ic_mode) == 'profile') THEN 1963 2282 output_var_table(21) % grid => scalar_profile_grid 1964 2283 output_var_table(21) % intermediate_grid => scalar_profile_intermediate … … 2026 2345 2027 2346 output_var_table(27) = init_nc_var( & 2028 name = 'init_ w',&2347 name = 'init_atmosphere_w', & 2029 2348 std_name = "", & 2030 2349 long_name = "initial wind component in z direction", & … … 2035 2354 grid = w_initial_grid, & 2036 2355 intermediate_grid = w_initial_intermediate, & 2037 is_profile = (TRIM( mode) == 'profile')&2038 ) 2039 IF (TRIM( mode) == 'profile') THEN2356 is_profile = (TRIM(ic_mode) == 'profile') & 2357 ) 2358 IF (TRIM(ic_mode) == 'profile') THEN 2040 2359 output_var_table(27) % grid => w_profile_grid 2041 2360 output_var_table(27) % intermediate_grid => w_profile_intermediate … … 2209 2528 intermediate_grid = palm_intermediate & 2210 2529 ) 2211 2530 ! 2531 !------------------------------------------------------------------------------ 2532 ! Section 2.2: Idealized large-scale forcings 2533 !------------------------------------------------------------------------------ 2212 2534 output_var_table(42) = init_nc_var( & 2213 2535 name = 'surface_forcing_surface_pressure', & … … 2227 2549 long_name = "geostrophic wind (u component)", & 2228 2550 units = "m/s", & 2229 kind = " profile",&2230 input_id = 1, & 2231 output_file = output_file, & 2232 grid = palm_grid,&2233 intermediate_grid = palm_intermediate&2551 kind = "constant scalar profile", & 2552 input_id = 1, & 2553 output_file = output_file, & 2554 grid = scalar_profile_grid, & 2555 intermediate_grid = scalar_profile_intermediate & 2234 2556 ) 2235 2557 … … 2239 2561 long_name = "geostrophic wind (v component)", & 2240 2562 units = "m/s", & 2241 kind = "profile", & 2242 input_id = 1, & 2243 output_file = output_file, & 2244 grid = palm_grid, & 2245 intermediate_grid = palm_intermediate & 2246 ) 2563 kind = "constant scalar profile", & 2564 input_id = 1, & 2565 output_file = output_file, & 2566 grid = scalar_profile_grid, & 2567 intermediate_grid = scalar_profile_intermediate & 2568 ) 2569 2570 output_var_table(45) = init_nc_var( & 2571 name = 'nudging_u', & 2572 std_name = "", & 2573 long_name = "wind component in x direction", & 2574 units = "m/s", & 2575 kind = "large-scale scalar forcing", & 2576 input_id = 1, & 2577 output_file = output_file, & 2578 grid = scalar_profile_grid, & 2579 intermediate_grid = scalar_profile_intermediate & 2580 ) 2581 output_var_table(45) % to_be_processed = ls_forcing_variables_required 2582 2583 output_var_table(46) = init_nc_var( & 2584 name = 'nudging_v', & 2585 std_name = "", & 2586 long_name = "wind component in y direction", & 2587 units = "m/s", & 2588 kind = "large-scale scalar forcing", & 2589 input_id = 1, & 2590 output_file = output_file, & 2591 grid = scalar_profile_grid, & 2592 intermediate_grid = scalar_profile_intermediate & 2593 ) 2594 output_var_table(46) % to_be_processed = ls_forcing_variables_required 2595 2596 output_var_table(47) = init_nc_var( & 2597 name = 'ls_forcing_sub_w', & 2598 std_name = "", & 2599 long_name = "subsidence velocity of w", & 2600 units = "m/s", & 2601 kind = "large-scale w forcing", & 2602 input_id = 1, & 2603 output_file = output_file, & 2604 grid = w_profile_grid, & 2605 intermediate_grid = w_profile_intermediate & 2606 ) 2607 output_var_table(47) % to_be_processed = ls_forcing_variables_required 2608 2609 output_var_table(48) = init_nc_var( & 2610 name = 'nudging_w', & 2611 std_name = "", & 2612 long_name = "wind component in w direction", & 2613 units = "m/s", & 2614 kind = "large-scale w forcing", & 2615 input_id = 1, & 2616 output_file = output_file, & 2617 grid = w_profile_grid, & 2618 intermediate_grid = w_profile_intermediate & 2619 ) 2620 output_var_table(48) % to_be_processed = ls_forcing_variables_required 2621 2622 2623 output_var_table(49) = init_nc_var( & 2624 name = 'ls_forcing_adv_pt', & 2625 std_name = "", & 2626 long_name = "advection of potential temperature", & 2627 units = "K/s", & 2628 kind = "large-scale scalar forcing", & 2629 input_id = 1, & 2630 output_file = output_file, & 2631 grid = scalar_profile_grid, & 2632 intermediate_grid = scalar_profile_intermediate & 2633 ) 2634 output_var_table(49) % to_be_processed = ls_forcing_variables_required 2635 2636 output_var_table(50) = init_nc_var( & 2637 name = 'ls_forcing_sub_pt', & 2638 std_name = "", & 2639 long_name = "subsidence velocity of potential temperature", & 2640 units = "K/s", & 2641 kind = "large-scale scalar forcing", & 2642 input_id = 1, & 2643 output_file = output_file, & 2644 grid = scalar_profile_grid, & 2645 intermediate_grid = scalar_profile_intermediate & 2646 ) 2647 output_var_table(50) % to_be_processed = ls_forcing_variables_required 2648 2649 output_var_table(51) = init_nc_var( & 2650 name = 'nudging_pt', & 2651 std_name = "", & 2652 long_name = "potential temperature", & 2653 units = "K", & 2654 kind = "large-scale scalar forcing", & 2655 input_id = 1, & 2656 output_file = output_file, & 2657 grid = scalar_profile_grid, & 2658 intermediate_grid = scalar_profile_intermediate & 2659 ) 2660 output_var_table(51) % to_be_processed = ls_forcing_variables_required 2661 2662 output_var_table(52) = init_nc_var( & 2663 name = 'ls_forcing_adv_qv', & 2664 std_name = "", & 2665 long_name = "advection of specific humidity", & 2666 units = "kg/kg/s", & 2667 kind = "large-scale scalar forcing", & 2668 input_id = 1, & 2669 output_file = output_file, & 2670 grid = scalar_profile_grid, & 2671 intermediate_grid = scalar_profile_intermediate & 2672 ) 2673 output_var_table(52) % to_be_processed = ls_forcing_variables_required 2674 2675 2676 output_var_table(53) = init_nc_var( & 2677 name = 'ls_forcing_sub_qv', & 2678 std_name = "", & 2679 long_name = "subsidence velocity of specific humidity", & 2680 units = "kg/kg/s", & 2681 kind = "large-scale scalar forcing", & 2682 input_id = 1, & 2683 output_file = output_file, & 2684 grid = scalar_profile_grid, & 2685 intermediate_grid = scalar_profile_intermediate & 2686 ) 2687 output_var_table(53) % to_be_processed = ls_forcing_variables_required 2688 2689 output_var_table(54) = init_nc_var( & 2690 name = 'nudging_qv', & 2691 std_name = "", & 2692 long_name = "specific humidity", & 2693 units = "kg/kg", & 2694 kind = "large-scale scalar forcing", & 2695 input_id = 1, & 2696 output_file = output_file, & 2697 grid = scalar_profile_grid, & 2698 intermediate_grid = scalar_profile_intermediate & 2699 ) 2700 output_var_table(54) % to_be_processed = ls_forcing_variables_required 2701 2702 output_var_table(55) = init_nc_var( & 2703 name = 'nudging_tau', & 2704 std_name = "", & 2705 long_name = "nudging relaxation time scale", & 2706 units = "s", & 2707 kind = "constant scalar profile", & 2708 input_id = 1, & 2709 output_file = output_file, & 2710 grid = scalar_profile_grid, & 2711 intermediate_grid = scalar_profile_intermediate & 2712 ) 2713 output_var_table(55) % to_be_processed = ls_forcing_variables_required 2714 2247 2715 2248 2716 ! Attributes shared among all variables 2249 2717 output_var_table(:) % source = nc_source_text 2718 2250 2719 2251 2720 END SUBROUTINE setup_variable_tables … … 2260 2729 !------------------------------------------------------------------------------! 2261 2730 FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, & 2262 grid, intermediate_grid, output_file, is_profile) RESULT(var) 2731 grid, intermediate_grid, output_file, is_profile & 2732 ) RESULT(var) 2263 2733 2264 2734 CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind … … 2367 2837 var % task = "average profile" 2368 2838 2369 CASE 2370 var % lod = 12839 CASE( 'surface forcing' ) 2840 var % lod = -1 2371 2841 var % ndim = 3 2372 2842 var % dimids(3) = output_file % dimid_time … … 2377 2847 var % task = "interpolate_2d" 2378 2848 2379 CASE 2380 var % lod = 22849 CASE( 'left scalar', 'right scalar') ! same as right 2850 var % lod = -1 2381 2851 var % ndim = 3 2382 2852 var % dimids(3) = output_file % dimid_time … … 2389 2859 var % task = "interpolate_3d" 2390 2860 2391 CASE 2392 var % lod = 22861 CASE( 'north scalar', 'south scalar') ! same as south 2862 var % lod = -1 2393 2863 var % ndim = 3 2394 2864 var % dimids(3) = output_file % dimid_time … … 2401 2871 var % task = "interpolate_3d" 2402 2872 2403 CASE 2404 var % lod = 22873 CASE( 'top scalar', 'top w' ) 2874 var % lod = -1 2405 2875 var % ndim = 3 2406 2876 var % dimids(3) = output_file % dimid_time … … 2413 2883 var % task = "interpolate_3d" 2414 2884 2415 CASE 2416 var % lod = 22885 CASE( 'left u', 'right u' ) 2886 var % lod = -1 2417 2887 var % ndim = 3 2418 2888 var % dimids(3) = output_file % dimid_time … … 2425 2895 var % task = "interpolate_3d" 2426 2896 2427 CASE 2428 var % lod = 22897 CASE( 'north u', 'south u' ) 2898 var % lod = -1 2429 2899 var % ndim = 3 2430 2900 var % dimids(3) = output_file % dimid_time !t … … 2437 2907 var % task = "interpolate_3d" 2438 2908 2439 CASE 2440 var % lod = 22909 CASE( 'top u' ) 2910 var % lod = -1 2441 2911 var % ndim = 3 2442 2912 var % dimids(3) = output_file % dimid_time !t … … 2449 2919 var % task = "interpolate_3d" 2450 2920 2451 CASE 2452 var % lod = 22921 CASE( 'left v', 'right v' ) 2922 var % lod = -1 2453 2923 var % ndim = 3 2454 2924 var % dimids(3) = output_file % dimid_time … … 2461 2931 var % task = "interpolate_3d" 2462 2932 2463 CASE 2464 var % lod = 22933 CASE( 'north v', 'south v' ) 2934 var % lod = -1 2465 2935 var % ndim = 3 2466 2936 var % dimids(3) = output_file % dimid_time !t … … 2473 2943 var % task = "interpolate_3d" 2474 2944 2475 CASE 2476 var % lod = 22945 CASE( 'top v' ) 2946 var % lod = -1 2477 2947 var % ndim = 3 2478 2948 var % dimids(3) = output_file % dimid_time !t … … 2485 2955 var % task = "interpolate_3d" 2486 2956 2487 CASE 2488 var % lod = 22957 CASE( 'left w', 'right w') 2958 var % lod = -1 2489 2959 var % ndim = 3 2490 2960 var % dimids(3) = output_file % dimid_time … … 2497 2967 var % task = "interpolate_3d" 2498 2968 2499 CASE 2500 var % lod = 22969 CASE( 'north w', 'south w' ) 2970 var % lod = -1 2501 2971 var % ndim = 3 2502 2972 var % dimids(3) = output_file % dimid_time !t … … 2509 2979 var % task = "interpolate_3d" 2510 2980 2511 CASE 2981 CASE( 'time series' ) 2512 2982 var % lod = 0 2513 2983 var % ndim = 1 … … 2517 2987 var % task = "average scalar" 2518 2988 2519 CASE ( 'profile' )2520 var % lod = 22989 CASE( 'constant scalar profile' ) 2990 var % lod = -1 2521 2991 var % ndim = 2 2522 2992 var % dimids(2) = output_file % dimid_time !t … … 2525 2995 var % dimvarids(1) = output_file % dimvarids_scl(3) 2526 2996 var % to_be_processed = .TRUE. 2527 var % task = "profile" 2997 var % task = "set profile" 2998 2999 CASE( 'large-scale scalar forcing' ) 3000 var % lod = -1 3001 var % ndim = 2 3002 var % dimids(2) = output_file % dimid_time !t 3003 var % dimids(1) = output_file % dimids_scl(3) !z 3004 var % dimvarids(2) = output_file % dimvarid_time 3005 var % dimvarids(1) = output_file % dimvarids_scl(3) 3006 var % to_be_processed = ls_forcing_variables_required 3007 var % task = "average large-scale profile" 3008 3009 CASE( 'large-scale w forcing' ) 3010 var % lod = -1 3011 var % ndim = 2 3012 var % dimids(2) = output_file % dimid_time !t 3013 var % dimids(1) = output_file % dimids_vel(3) !z 3014 var % dimvarids(2) = output_file % dimvarid_time 3015 var % dimvarids(1) = output_file % dimvarids_vel(3) 3016 var % to_be_processed = ls_forcing_variables_required 3017 var % task = "average large-scale profile" 2528 3018 2529 3019 CASE DEFAULT … … 2560 3050 2561 3051 2562 SUBROUTINE input_file_list(start_date_string, start_hour, end_hour, &2563 step_hour, path, prefix, suffix, file_list)3052 SUBROUTINE get_input_file_list(start_date_string, start_hour, end_hour, & 3053 step_hour, path, prefix, suffix, file_list) 2564 3054 2565 3055 CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string … … 2568 3058 CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:) 2569 3059 2570 INTEGER :: number_of_ files, hour, i3060 INTEGER :: number_of_intervals, hour, i 2571 3061 CHARACTER(LEN=DATE) :: date_string 2572 3062 2573 number_of_files = end_hour - start_hour + 1 2574 2575 ALLOCATE( file_list(number_of_files) ) 2576 2577 DO i = 1, number_of_files 2578 hour = start_hour + (i-1) * step_hour 3063 number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour ) 3064 ALLOCATE( file_list(number_of_intervals + 1) ) 3065 3066 DO i = 0, number_of_intervals 3067 hour = start_hour + i * step_hour 2579 3068 date_string = add_hours_to(start_date_string, hour) 2580 3069 2581 file_list(i ) = TRIM(path) // TRIM(prefix) // TRIM(date_string) // &2582 TRIM(suffix) // '.nc'2583 message = "Set up input file name '" // TRIM(file_list(i )) //"'"3070 file_list(i+1) = TRIM(path) // TRIM(prefix) // TRIM(date_string) // & 3071 TRIM(suffix) // '.nc' 3072 message = "Set up input file name '" // TRIM(file_list(i+1)) //"'" 2584 3073 CALL report('input_file_list', message) 2585 3074 END DO 2586 3075 2587 END SUBROUTINE input_file_list3076 END SUBROUTINE get_input_file_list 2588 3077 2589 3078 … … 2607 3096 2608 3097 REAL(dp), ALLOCATABLE :: basic_state_pressure(:) 2609 TYPE(container), ALLOCATABLE :: compute_buffer(:)3098 TYPE(container), ALLOCATABLE :: preprocess_buffer(:) 2610 3099 INTEGER :: hour, dt 2611 3100 INTEGER :: i, j, k … … 2617 3106 2618 3107 CASE( 'velocities' ) 2619 ! Allocate a compute puffer with the same number of arrays as the input2620 ALLOCATE( compute_buffer( SIZE(input_buffer) ) )3108 ! Allocate a compute buffer with the same number of arrays as the input 3109 ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) ) 2621 3110 2622 3111 ! Allocate u and v arrays with scalar dimensions … … 2624 3113 ny = SIZE(input_buffer(1) % array, 2) 2625 3114 nz = SIZE(input_buffer(1) % array, 3) 2626 ALLOCATE( compute_buffer(1) % array(nx, ny, nz) ) ! u buffer2627 ALLOCATE( compute_buffer(2) % array(nx, ny, nz) ) ! v buffer3115 ALLOCATE( preprocess_buffer(1) % array(nx, ny, nz) ) ! u buffer 3116 ALLOCATE( preprocess_buffer(2) % array(nx, ny, nz) ) ! v buffer 2628 3117 2629 3118 CALL run_control('time', 'alloc') … … 2632 3121 CALL centre_velocities( u_face = input_buffer(1) % array, & 2633 3122 v_face = input_buffer(2) % array, & 2634 u_centre = compute_buffer(1) % array,&2635 v_centre = compute_buffer(2) % array )3123 u_centre = preprocess_buffer(1) % array, & 3124 v_centre = preprocess_buffer(2) % array ) 2636 3125 2637 ! rotate U and V to PALM-4U orientation and overwrite U and V with 2638 ! rotated velocities 2639 DO k = 1, nz 2640 DO j = 2, ny 2641 DO i = 2, nx 2642 CALL uv2uvrot( urot = compute_buffer(1) % array(i,j,k), & 2643 vrot = compute_buffer(2) % array(i,j,k), & 2644 rlat = cosmo_grid % lat(j-1), & 2645 rlon = cosmo_grid % lon(i-1), & 2646 pollat = phi_cn, & 2647 pollon = lambda_cn, & 2648 u = input_buffer(1) % array(i,j,k), & 2649 v = input_buffer(2) % array(i,j,k) ) 2650 END DO 2651 END DO 2652 END DO 3126 cfg % rotation_method = 'rotated-pole' 3127 SELECT CASE(cfg % rotation_method) 3128 3129 CASE('rotated-pole') 3130 ! rotate U and V to PALM-4U orientation and overwrite U and V with 3131 ! rotated velocities 3132 DO k = 1, nz 3133 DO j = 2, ny 3134 DO i = 2, nx 3135 CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k), & 3136 vrot = preprocess_buffer(2) % array(i,j,k), & 3137 rlat = cosmo_grid % lat(j-1), & 3138 rlon = cosmo_grid % lon(i-1), & 3139 pollat = phi_cn, & 3140 pollon = lambda_cn, & 3141 u = input_buffer(1) % array(i,j,k), & 3142 v = input_buffer(2) % array(i,j,k) ) 3143 END DO 3144 END DO 3145 END DO 3146 3147 CASE DEFAULT 3148 message = "Rotation method '" // TRIM(cfg % rotation_method) // & 3149 "' not recognized." 3150 CALL abort('preprocess', message) 3151 3152 END SELECT 3153 3154 ! set values 2653 3155 input_buffer(1) % array(1,:,:) = 0.0_dp 2654 3156 input_buffer(2) % array(1,:,:) = 0.0_dp … … 2659 3161 CALL run_control('time', 'comp') 2660 3162 2661 DEALLOCATE( compute_buffer )3163 DEALLOCATE( preprocess_buffer ) 2662 3164 CALL run_control('time', 'alloc') 2663 3165
Note: See TracChangeset
for help on using the changeset viewer.
