Changeset 3182 for palm/trunk/UTIL/inifor/src
- Timestamp:
- Jul 27, 2018 1:36:03 PM (6 years ago)
- Location:
- palm/trunk/UTIL/inifor/src
- Files:
-
- 1 added
- 8 edited
-
control.f90 (modified) (3 diffs)
-
defs.f90 (modified) (3 diffs)
-
get_netcdf_variable.inc (added)
-
grid.f90 (modified) (93 diffs)
-
inifor.f90 (modified) (14 diffs)
-
io.f90 (modified) (21 diffs)
-
transform.f90 (modified) (16 diffs)
-
types.f90 (modified) (8 diffs)
-
util.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/control.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Added version and copyright output 23 24 ! 24 25 ! … … 42 43 43 44 USE defs, & 44 ONLY: LNAME, dp 45 ONLY: LNAME, dp, VERSION, COPYRIGHT 45 46 46 47 USE util, & … … 95 96 96 97 END SUBROUTINE abort 98 99 100 SUBROUTINE print_version() 101 PRINT *, "INIFOR " // VERSION 102 PRINT *, COPYRIGHT 103 END SUBROUTINE print_version 97 104 98 105 -
palm/trunk/UTIL/inifor/src/defs.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Updated defaults for soil extrapolation steps and nudging time-scale 24 ! Improved handling of the start date string 25 ! Added gas constant for water vapor 26 ! Bumped INIFOR version number 23 27 ! 24 28 ! … … 62 66 REAL(dp), PARAMETER :: T_SL = 288.15_dp !< Reference temperature for computation of COSMO-DE's basic state pressure [K] 63 67 REAL(dp), PARAMETER :: BETA = 42.0_dp !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic state pressure [K] 64 REAL(dp), PARAMETER :: RD = 287.05_dp !< specific gar constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K] 68 REAL(dp), PARAMETER :: RD = 287.05_dp !< specific gas constant of dry air, used in computation of COSMO-DE's basic state [J/kg/K] 69 REAL(dp), PARAMETER :: RV = 461.51_dp !< specific gas constant of water vapor [J/kg/K] 65 70 REAL(dp), PARAMETER :: G = 9.80665_dp !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s] 66 71 REAL(dp), PARAMETER :: RHO_L = 1e3_dp !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3] … … 72 77 73 78 ! INIFOR parameters 74 INTEGER, PARAMETER :: FILL_ITERATIONS = 10 !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-] 75 REAL(dp), PARAMETER :: FORCING_FREQ = 3600.0_dp !< Reference pressure for potential temperature [Pa] 76 CHARACTER(LEN=*), PARAMETER :: VERSION = '1.1.4' !< path to script for generating input file names 79 INTEGER, PARAMETER :: FILL_ITERATIONS = 5 !< Number of iterations for extrapolating soil data into COSMO-DE water cells [-] 80 INTEGER, PARAMETER :: FORCING_STEP = 1 !< Number of hours between forcing time steps [h] 81 REAL(dp), PARAMETER :: NUDGING_TAU = 21600.0_dp !< Nudging relaxation time scale [s] 82 CHARACTER(LEN=*), PARAMETER :: VERSION = '1.3.0' !< INIFOR version number 83 CHARACTER(LEN=*), PARAMETER :: COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // & 84 NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice 77 85 78 86 END MODULE defs -
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 -
palm/trunk/UTIL/inifor/src/inifor.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Introduced new PALM grid stretching 24 ! Renamend initial-condition mode variable 'mode' to 'ic_mode' 25 ! Improved log messages 23 26 ! 24 27 ! … … 48 51 ONLY: setup_parameters, setup_grids, setup_variable_tables, & 49 52 setup_io_groups, fini_grids, fini_variables, fini_io_groups, & 50 fini_file_lists, preprocess, 53 fini_file_lists, preprocess, origin_lon, origin_lat, & 51 54 output_file, io_group_list, output_var_table, & 52 cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, mode,&53 imin, imax, jmin,jmax55 cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg, & 56 average_imin, average_imax, average_jmin, average_jmax 54 57 55 58 USE io … … 84 87 85 88 ! Initialize the netCDF output file and define dimensions 86 CALL setup_netcdf_dimensions(output_file, palm_grid) 89 CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date, & 90 origin_lon, origin_lat) 87 91 CALL run_control('time', 'write') 88 92 89 93 ! Set up the tables containing the input and output variables and set 90 94 ! the corresponding netCDF dimensions for each output variable 91 CALL setup_variable_tables( mode)95 CALL setup_variable_tables(cfg % ic_mode) 92 96 CALL run_control('time', 'write') 93 97 … … 95 99 CALL setup_netcdf_variables(output_file % name, output_var_table) 96 100 97 CALL setup_io_groups() 101 CALL setup_io_groups() 98 102 CALL run_control('time', 'init') 99 103 … … 118 122 CALL run_control('time', 'comp') 119 123 124 !TODO: move this assertion into 'preprocess'. 120 125 IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) ) THEN 121 126 message = "Input buffers for group '" // TRIM(group % kind) // & … … 159 164 CASE DEFAULT 160 165 161 CALL abort("main loop", 'Not a soil variable') 166 message = "'" // TRIM(output_var % kind) // "' is not a soil variable" 167 CALL abort("main loop", message) 162 168 163 169 END SELECT … … 173 179 ALLOCATE( output_arr( 0:output_var % grid % nx, & 174 180 0:output_var % grid % ny, & 175 0:output_var % grid % nz ) )181 1:output_var % grid % nz ) ) 176 182 177 183 CALL run_control('time', 'alloc') … … 187 193 ALLOCATE( output_arr( 0:output_var % grid % nx, & 188 194 0:output_var % grid % ny, & 189 0:output_var % grid % nz ) )195 1:output_var % grid % nz ) ) 190 196 CALL run_control('time', 'alloc') 191 197 … … 193 199 CALL average_profile( & 194 200 input_buffer(output_var % input_id) % array(:,:,:), & 195 output_arr(:,:,:), imin, imax, jmin, jmax, & 201 output_arr(:,:,:), average_imin, average_imax, & 202 average_jmin, average_jmax, & 196 203 output_var % intermediate_grid, & 197 204 output_var % grid) … … 205 212 CALL run_control('time', 'comp') 206 213 207 CASE ( ' profile' )214 CASE ( 'set profile' ) 208 215 209 ALLOCATE( output_arr( 1, 1, 0:nz ) )216 ALLOCATE( output_arr( 1, 1, 1:nz ) ) 210 217 CALL run_control('time', 'alloc') 211 218 … … 217 224 CASE('ls_forcing_vg') 218 225 output_arr(1, 1, :) = vg 226 227 CASE('nudging_tau') 228 output_arr(1, 1, :) = NUDGING_TAU 219 229 220 230 CASE DEFAULT … … 225 235 END SELECT 226 236 CALL run_control('time', 'comp') 237 238 CASE('average large-scale profile') 239 message = "Averaging of large-scale forcing profiles " //& 240 "has not been implemented, yet." 241 CALL abort('main loop', message) 242 !ALLOCATE( output_arr( 1, 1, 1:nz ) ) 227 243 228 244 CASE DEFAULT … … 269 285 ELSE 270 286 271 message = "Skipping IO group '" // TRIM(group % kind) // "'"287 message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'" 272 288 IF ( ALLOCATED(group % in_var_list) ) THEN 273 289 message = TRIM(message) // " with input variable '" // & … … 291 307 CALL run_control('report', 'void') 292 308 293 message = "Finished writing forcing file'" // TRIM(output_file % name) // &309 message = "Finished writing dynamic driver '" // TRIM(output_file % name) // & 294 310 "' successfully." 295 311 CALL report('main loop', message) -
palm/trunk/UTIL/inifor/src/io.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Introduced new PALM grid stretching 24 ! Updated variable names and metadata for PIDS v1.9 compatibility 25 ! Improved handling of the start date string 26 ! Better compatibility with older Intel compilers: 27 ! - avoiding implicit array allocation with new get_netcdf_variable() 28 ! subroutine instead of function 29 ! Improved command line interface: 30 ! - Added configuration validation 31 ! - New options to configure input file prefixes 32 ! - GNU-style short and long option names 33 ! - Added version and copyright output 23 34 ! 24 ! 35 ! 25 36 ! Former revisions: 26 37 ! ----------------- … … 43 54 USE control 44 55 USE defs, & 45 ONLY: DATE, SNAME, PATH, PI, dp, TO_RADIANS, TO_DEGREES, VERSION56 ONLY: DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION 46 57 USE netcdf 47 58 USE types 48 59 USE util, & 49 ONLY: reverse, str 60 ONLY: reverse, str, real_to_str 50 61 51 62 IMPLICIT NONE 52 63 64 INTERFACE get_netcdf_variable 65 MODULE PROCEDURE get_netcdf_variable_int 66 MODULE PROCEDURE get_netcdf_variable_real 67 END INTERFACE get_netcdf_variable 68 69 PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real 70 53 71 CONTAINS 72 73 SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer) 74 75 CHARACTER(LEN=PATH), INTENT(IN) :: in_file 76 TYPE(nc_var), INTENT(INOUT) :: in_var 77 INTEGER(hp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:) 78 79 INCLUDE 'get_netcdf_variable.inc' 80 81 END SUBROUTINE get_netcdf_variable_int 82 83 84 SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer) 85 86 CHARACTER(LEN=PATH), INTENT(IN) :: in_file 87 TYPE(nc_var), INTENT(INOUT) :: in_var 88 REAL(dp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:) 89 90 INCLUDE 'get_netcdf_variable.inc' 91 92 END SUBROUTINE get_netcdf_variable_real 93 54 94 55 95 SUBROUTINE netcdf_define_variable(var, ncid) … … 59 99 60 100 CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT, var % dimids(1:var % ndim), var % varid)) 61 CALL check(nf90_put_att(ncid, var % varid, "standard_name", var % standard_name))62 101 CALL check(nf90_put_att(ncid, var % varid, "long_name", var % long_name)) 63 102 CALL check(nf90_put_att(ncid, var % varid, "units", var % units)) 64 CALL check(nf90_put_att(ncid, var % varid, "lod", var % lod)) 103 IF ( var % lod .GE. 0 ) THEN 104 CALL check(nf90_put_att(ncid, var % varid, "lod", var % lod)) 105 END IF 65 106 CALL check(nf90_put_att(ncid, var % varid, "source", var % source)) 66 107 CALL check(nf90_put_att(ncid, var % varid, "_FillValue", NF90_FILL_REAL)) … … 94 135 !> parameters for the PALM-4U computational grid. 95 136 !------------------------------------------------------------------------------! 96 SUBROUTINE parse_command_line_arguments( start_date, hhl_file, & 97 soiltyp_file, static_driver_file, input_path, output_file, & 98 namelist_file, ug, vg, p0, z0, mode ) 99 100 CHARACTER(LEN=PATH), INTENT(INOUT) :: hhl_file, soiltyp_file, & 101 static_driver_file, input_path, output_file, namelist_file 102 CHARACTER(LEN=SNAME), INTENT(INOUT) :: mode 103 REAL(dp), INTENT(INOUT) :: ug, vg, p0, z0 104 CHARACTER(LEN=DATE), INTENT(INOUT) :: start_date 105 106 CHARACTER(LEN=PATH) :: option, arg 107 INTEGER :: arg_count, i 137 SUBROUTINE parse_command_line_arguments( cfg ) 138 139 TYPE(inifor_config), INTENT(INOUT) :: cfg 140 141 CHARACTER(LEN=PATH) :: option, arg 142 INTEGER :: arg_count, i 108 143 109 144 arg_count = COMMAND_ARGUMENT_COUNT() … … 111 146 112 147 ! Every option should have an argument. 113 IF ( MOD(arg_count, 2) .NE. 0 ) THEN114 message = "Syntax error in command line."115 CALL abort('parse_command_line_arguments', message)116 END IF148 !IF ( MOD(arg_count, 2) .NE. 0 ) THEN 149 ! message = "Syntax error in command line." 150 ! CALL abort('parse_command_line_arguments', message) 151 !END IF 117 152 118 153 message = "The -clon and -clat command line options are depricated. " // & 119 154 "Please remove them form your inifor command and specify the " // & 120 155 "location of the PALM-4U origin either" // NEW_LINE(' ') // & 121 " - by setting the namelist parameters ' origin_lon' and 'origin_lat, or'" // NEW_LINE(' ') // &156 " - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // & 122 157 " - by providing a static driver netCDF file via the -static command-line option." 123 158 124 ! Loop through option/argument pairs.125 DO i = 1, arg_count, 2159 i = 1 160 DO WHILE (i .LE. arg_count) 126 161 127 162 CALL GET_COMMAND_ARGUMENT( i, option ) 128 CALL GET_COMMAND_ARGUMENT( i+1, arg )129 163 130 164 SELECT CASE( TRIM(option) ) 131 165 132 CASE( '-date' ) 133 start_date = TRIM(arg) 166 CASE( '-date', '-d', '--date' ) 167 CALL get_option_argument( i, arg ) 168 cfg % start_date = TRIM(arg) 134 169 135 170 ! Elevation of the PALM-4U domain above sea level 136 CASE( '-z0' ) 137 READ(arg, *) z0 171 CASE( '-z0', '-z', '--elevation' ) 172 CALL get_option_argument( i, arg ) 173 READ(arg, *) cfg % z0 138 174 139 175 ! surface pressure, at z0 140 CASE( '-p0' ) 141 READ(arg, *) p0 142 143 ! surface pressure, at z0 144 CASE( '-ug' ) 145 READ(arg, *) ug 146 147 ! surface pressure, at z0 148 CASE( '-vg' ) 149 READ(arg, *) vg 150 151 ! Domain centre geographical longitude 152 CASE( '-clon' ) 176 CASE( '-p0', '-r', '--surface-pressure' ) 177 CALL get_option_argument( i, arg ) 178 READ(arg, *) cfg % p0 179 180 ! geostrophic wind in x direction 181 CASE( '-ug', '-u', '--geostrophic-u' ) 182 CALL get_option_argument( i, arg ) 183 READ(arg, *) cfg % ug 184 185 ! geostrophic wind in y direction 186 CASE( '-vg', '-v', '--geostrophic-v' ) 187 CALL get_option_argument( i, arg ) 188 READ(arg, *) cfg % vg 189 190 ! domain centre geographical longitude and latitude 191 CASE( '-clon', '-clat' ) 153 192 CALL abort('parse_command_line_arguments', message) 154 193 !READ(arg, *) lambda_cg 155 194 !lambda_cg = lambda_cg * TO_RADIANS 156 157 ! Domain centre geographical latitude158 CASE( '-clat' )159 CALL abort('parse_command_line_arguments', message)160 195 !READ(arg, *) phi_cg 161 196 !phi_cg = phi_cg * TO_RADIANS 162 197 163 CASE( '-path' ) 164 input_path = TRIM(arg) 165 166 CASE( '-hhl' ) 167 hhl_file = TRIM(arg) 168 169 CASE( '-static' ) 170 static_driver_file = TRIM(arg) 171 172 CASE( '-soil' ) 173 soiltyp_file = TRIM(arg) 174 175 CASE( '-o' ) 176 output_file = TRIM(arg) 177 178 CASE( '-n' ) 179 namelist_file = TRIM(arg) 180 181 ! Initialization mode: 'profile' / 'volume' 182 CASE( '-mode' ) 183 mode = TRIM(arg) 184 185 SELECT CASE( TRIM(mode) ) 186 187 CASE( 'profile' ) 188 189 CASE DEFAULT 190 message = "Mode '" // TRIM(mode) // "' is not supported. " //& 191 "Currently, '-mode profile' is the only supported option. " //& 192 "Select this one or omit the -mode option entirely." 193 CALL abort( 'parse_command_line_arguments', message ) 194 END SELECT 198 CASE( '-path', '-p', '--path' ) 199 CALL get_option_argument( i, arg ) 200 cfg % input_path = TRIM(arg) 201 202 CASE( '-hhl', '-l', '--hhl-file' ) 203 CALL get_option_argument( i, arg ) 204 cfg % hhl_file = TRIM(arg) 205 206 CASE( '-static', '-t', '--static-driver' ) 207 CALL get_option_argument( i, arg ) 208 cfg % static_driver_file = TRIM(arg) 209 210 CASE( '-soil', '-s', '--soil-file') 211 CALL get_option_argument( i, arg ) 212 cfg % soiltyp_file = TRIM(arg) 213 214 CASE( '--flow-prefix') 215 CALL get_option_argument( i, arg ) 216 cfg % flow_prefix = TRIM(arg) 217 218 CASE( '--radiation-prefix') 219 CALL get_option_argument( i, arg ) 220 cfg % radiation_prefix = TRIM(arg) 221 222 CASE( '--soil-prefix') 223 CALL get_option_argument( i, arg ) 224 cfg % soil_prefix = TRIM(arg) 225 226 CASE( '--soilmoisture-prefix') 227 CALL get_option_argument( i, arg ) 228 cfg % soilmoisture_prefix = TRIM(arg) 229 230 CASE( '-o', '--output' ) 231 CALL get_option_argument( i, arg ) 232 cfg % output_file = TRIM(arg) 233 234 CASE( '-n', '--namelist' ) 235 CALL get_option_argument( i, arg ) 236 cfg % namelist_file = TRIM(arg) 237 238 ! initial condition mode: 'profile' / 'volume' 239 CASE( '-mode', '-i', '--init-mode' ) 240 CALL get_option_argument( i, arg ) 241 cfg % ic_mode = TRIM(arg) 242 243 ! boundary conditions / forcing mode: 'ideal' / 'real' 244 CASE( '-f', '--forcing-mode' ) 245 CALL get_option_argument( i, arg ) 246 cfg % bc_mode = TRIM(arg) 247 248 CASE( '--version' ) 249 CALL print_version() 250 STOP 251 252 CASE( '--help' ) 253 CALL print_version() 254 PRINT *, "" 255 PRINT *, "For a list of command-line options have a look at the README file." 256 STOP 195 257 196 258 CASE DEFAULT 197 message = "unknown option '" // TRIM(option (2:)) // "'."259 message = "unknown option '" // TRIM(option) // "'." 198 260 CALL abort('parse_command_line_arguments', message) 199 261 200 262 END SELECT 263 264 i = i + 1 201 265 202 266 END DO … … 210 274 211 275 END SUBROUTINE parse_command_line_arguments 276 277 278 SUBROUTINE get_option_argument(i, arg) 279 CHARACTER(LEN=PATH), INTENT(INOUT) :: arg 280 INTEGER, INTENT(INOUT) :: i 281 282 i = i + 1 283 CALL GET_COMMAND_ARGUMENT(i, arg) 284 285 END SUBROUTINE 286 287 288 SUBROUTINE validate_config(cfg) 289 TYPE(inifor_config), INTENT(IN) :: cfg 290 LOGICAL :: all_files_present 291 292 all_files_present = .TRUE. 293 all_files_present = all_files_present .AND. file_present(cfg % hhl_file) 294 all_files_present = all_files_present .AND. file_present(cfg % namelist_file) 295 all_files_present = all_files_present .AND. file_present(cfg % output_file) 296 all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file) 297 298 ! Only check optional static driver file name, if it has been given. 299 IF (TRIM(cfg % static_driver_file) .NE. '') THEN 300 all_files_present = all_files_present .AND. file_present(cfg % static_driver_file) 301 END IF 302 303 IF (.NOT. all_files_present) THEN 304 message = "INIFOR configuration invalid; some input files are missing." 305 CALL abort( 'validate_config', message ) 306 END IF 307 308 309 SELECT CASE( TRIM(cfg % ic_mode) ) 310 CASE( 'profile', 'volume') 311 CASE DEFAULT 312 message = "Initialization mode '" // TRIM(cfg % ic_mode) //& 313 "' is not supported. " //& 314 "Please select either 'profile' or 'volume', " //& 315 "or omit the -i/--init-mode/-mode option entirely, which corresponds "//& 316 "to the latter." 317 CALL abort( 'validate_config', message ) 318 END SELECT 319 320 321 SELECT CASE( TRIM(cfg % bc_mode) ) 322 CASE( 'real', 'ideal') 323 CASE DEFAULT 324 message = "Forcing mode '" // TRIM(cfg % bc_mode) //& 325 "' is not supported. " //& 326 "Please select either 'real' or 'ideal', " //& 327 "or omit the -f/--forcing-mode option entirely, which corresponds "//& 328 "to the latter." 329 CALL abort( 'validate_config', message ) 330 END SELECT 331 332 333 END SUBROUTINE validate_config 334 335 336 LOGICAL FUNCTION file_present(filename) 337 CHARACTER(LEN=PATH), INTENT(IN) :: filename 338 339 INQUIRE(FILE=filename, EXIST=file_present) 340 341 IF (.NOT. file_present) THEN 342 message = "The given file '" // "' does not exist." 343 CALL report('file_present', message) 344 END IF 345 346 END FUNCTION file_present 212 347 213 348 … … 222 357 !> writes the actual data. 223 358 !------------------------------------------------------------------------------! 224 SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid) 359 SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid, & 360 start_date_string, origin_lon, origin_lat) 225 361 226 362 TYPE(nc_file), INTENT(INOUT) :: output_file 227 363 TYPE(grid_definition), INTENT(IN) :: palm_grid 228 229 CHARACTER (LEN=SNAME) :: date 364 CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string 365 REAL(dp), INTENT(IN) :: origin_lon, origin_lat 366 367 CHARACTER (LEN=8) :: date_string 368 CHARACTER (LEN=10) :: time_string 369 CHARACTER (LEN=5) :: zone_string 370 CHARACTER (LEN=SNAME) :: history_string 230 371 INTEGER :: ncid, nx, ny, nz, nt, dimids(3), dimvarids(3) 231 372 REAL(dp) :: z0 232 373 374 message = "Initializing PALM-4U dynamic driver file '" // & 375 TRIM(output_file % name) // "' and setting up dimensions." 376 CALL report('setup_netcdf_dimensions', message) 377 233 378 ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode. 379 #if defined( __netcdf4 ) 234 380 CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid)) 381 #else 382 CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid)) 383 #endif 235 384 236 385 ! … … 238 387 !- Section 1: Write global NetCDF attributes 239 388 !------------------------------------------------------------------------------ 240 CALL date_and_time(date) 389 CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string) 390 history_string = & 391 'Created on '// date_string // & 392 ' at ' // time_string(1:2) // ':' // time_string(3:4) // & 393 ' (UTC' // zone_string // ')' 394 241 395 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...')) 242 396 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach')) 243 397 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de')) 244 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', 'Created on '//date))398 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string))) 245 399 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--')) 246 400 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--')) 247 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', '--'))248 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', '--'))401 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)')))) 402 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)')))) 249 403 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION))) 250 404 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--')) … … 267 421 CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) ) 268 422 CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) ) 269 CALL check( nf90_def_dim(ncid, "z", nz +1, dimids(3)) )423 CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) ) 270 424 output_file % dimids_scl = dimids ! save dimids for later 271 425 … … 285 439 286 440 ! overwrite third dimid with the one of depth 287 CALL check(nf90_def_dim(ncid, " depth", SIZE(palm_grid % depths), dimids(3)) )441 CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) ) 288 442 output_file % dimids_soil = dimids ! save dimids for later 289 443 290 444 ! overwrite third dimvarid with the one of depth 291 CALL check(nf90_def_var(ncid, " depth", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))445 CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3))) 292 446 CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land")) 293 447 CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down")) … … 301 455 CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) ) 302 456 CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) ) 303 CALL check(nf90_def_dim(ncid, "zw", nz , dimids(3)) )457 CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) ) 304 458 output_file % dimids_vel = dimids ! save dimids for later 305 459 … … 328 482 CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time")) 329 483 CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time")) 330 CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", "seconds since...")) 484 CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", & 485 "seconds since " // start_date_string // " UTC")) 331 486 332 487 CALL check(nf90_enddef(ncid)) … … 363 518 INTEGER :: i, ncid 364 519 365 message = " Initializing PALM-4U forcing file'" // TRIM(filename) // "'."520 message = "Defining variables in dynamic driver '" // TRIM(filename) // "'." 366 521 CALL report('setup_netcdf_variables', message) 367 522 … … 374 529 375 530 IF ( var % to_be_processed ) THEN 376 message = " Definingvariable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."531 message = " variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'." 377 532 CALL report('setup_netcdf_variables', message) 378 533 … … 386 541 CALL check(nf90_close(ncid)) 387 542 388 message = " Forcing file'" // TRIM(filename) // "' initialized successfully."543 message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully." 389 544 CALL report('setup_netcdf_variables', message) 390 545 … … 447 602 448 603 input_var => group % in_var_list(1) 449 buffer(buf_id) % array = get_netcdf_variable( input_file, input_var )604 CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array) 450 605 CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name)) 451 606 … … 472 627 END IF 473 628 474 buffer(ivar) % array = get_netcdf_variable( input_file, input_var )629 CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array) 475 630 476 631 IF ( input_var % is_upside_down ) CALL reverse(buffer(ivar) % array) … … 545 700 546 701 CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value)) 702 CALL check(nf90_close(ncid)) 547 703 548 704 ELSE … … 555 711 556 712 END FUNCTION get_netcdf_attribute 557 558 559 560 FUNCTION get_netcdf_variable(in_file, in_var) RESULT(buffer)561 562 CHARACTER(LEN=PATH), INTENT(IN) :: in_file563 TYPE(nc_var), INTENT(INOUT) :: in_var564 REAL(dp), ALLOCATABLE :: buffer(:,:,:)565 INTEGER :: i, ncid, start(3)566 567 568 ! Read in_var NetCDF attributes569 IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &570 nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR ) THEN571 572 CALL check(nf90_get_att(ncid, in_var % varid, "long_name", in_var % long_name))573 CALL check(nf90_get_att(ncid, in_var % varid, "units", in_var % units))574 575 ! Read in_var NetCDF dimensions576 CALL check(nf90_inquire_variable( ncid, in_var % varid, &577 ndims = in_var % ndim, &578 dimids = in_var % dimids ))579 580 DO i = 1, in_var % ndim581 CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i), &582 name = in_var % dimname(i), &583 len = in_var % dimlen(i) ))584 END DO585 586 start = (/ 1, 1, 1 /)587 IF ( TRIM(in_var % name) .EQ. 'T_SO' ) THEN588 ! Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8589 in_var % dimlen(3) = in_var % dimlen(3) - 1590 591 ! Start reading from second level, e.g. depth = 0.005 instead of 0.0592 start(3) = 2593 END IF594 595 SELECT CASE(in_var % ndim)596 597 CASE (2)598 599 ALLOCATE( buffer( in_var % dimlen(1), &600 in_var % dimlen(2), &601 1 ) )602 603 CASE (3)604 605 ALLOCATE( buffer( in_var % dimlen(1), &606 in_var % dimlen(2), &607 in_var % dimlen(3) ) )608 CASE (4)609 610 ALLOCATE( buffer( in_var % dimlen(1), &611 in_var % dimlen(2), &612 in_var % dimlen(3) ) )613 CASE DEFAULT614 615 message = "Failed reading NetCDF variable " // &616 TRIM(in_var % name) // " with " // TRIM(str(in_var%ndim)) // &617 " dimensions because only two- and and three-dimensional" // &618 " variables are supported."619 CALL abort('get_netcdf_variable', message)620 621 END SELECT622 CALL run_control('time', 'alloc')623 624 ! TODO: Check for matching dimensions of buffer and var625 CALL check(nf90_get_var( ncid, in_var % varid, buffer, &626 start = start, &627 count = in_var % dimlen(1:3) ) )628 629 CALL run_control('time', 'read')630 ELSE631 632 message = "Failed to read '" // TRIM(in_var % name) // &633 "' from file '" // TRIM(in_file) // "'."634 CALL report('get_netcdf_variable', message)635 636 END IF637 638 CALL check(nf90_close(ncid))639 640 CALL run_control('time', 'read')641 642 END FUNCTION get_netcdf_variable643 713 644 714 … … 657 727 658 728 ! Skip time dimension for output 659 IF ( var_is_time_dependent ) THEN 660 ndim = var % ndim - 1 661 ELSE 662 ndim = var % ndim 663 END IF 729 ndim = var % ndim 730 IF ( var_is_time_dependent ) ndim = var % ndim - 1 664 731 665 732 start(:) = (/1,1,1,1/) … … 733 800 start=start(1:ndim+1) ) ) 734 801 735 CASE ( ' profile' )802 CASE ( 'constant scalar profile' ) 736 803 737 804 CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), & 738 805 start=start(1:ndim+1), & 739 806 count=count(1:ndim) ) ) 807 808 CASE ( 'large-scale scalar forcing', 'large-scale w forcing' ) 809 810 message = "Doing nothing in terms of writing large-scale forings." 811 CALL report('update_output', message) 740 812 741 813 CASE DEFAULT -
palm/trunk/UTIL/inifor/src/transform.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Introduced new PALM grid stretching 24 ! Removed unnecessary subroutine parameters 25 ! Renamed kcur to k_intermediate 23 26 ! 24 27 ! … … 80 83 TYPE(grid_definition), INTENT(IN) :: outgrid 81 84 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) 82 REAL(dp), INTENT(OUT) :: out_arr(0:,0:,0:) 83 84 INTEGER :: i, j, k, l, nx, ny, nz 85 86 nx = UBOUND(out_arr, 1) 87 ny = UBOUND(out_arr, 2) 85 REAL(dp), INTENT(OUT) :: out_arr(0:,0:,:) 86 87 INTEGER :: i, j, k, l, nz 88 88 89 nz = UBOUND(out_arr, 3) 89 90 90 DO j = 0, ny91 DO i = 0, nx92 DO k = nz, 0, -191 DO j = LBOUND(out_arr, 2), UBOUND(out_arr, 2) 92 DO i = LBOUND(out_arr, 1), UBOUND(out_arr, 1) 93 DO k = nz, LBOUND(out_arr, 3), -1 93 94 94 95 ! TODO: Remove IF clause and extrapolate based on a critical vertical … … 101 102 out_arr(i,j,k) = 0.0_dp 102 103 DO l = 1, 2 103 out_arr(i,j,k) = out_arr(i,j,k) + 104 outgrid % w_verti(i,j,k,l) * 104 out_arr(i,j,k) = out_arr(i,j,k) + & 105 outgrid % w_verti(i,j,k,l) * & 105 106 in_arr(i,j,outgrid % kk(i,j,k, l) ) 106 107 END DO … … 139 140 ! I index 0-based for the indices of the outvar to be consistent with the 140 141 ! outgrid indices and interpolation weights. 141 TYPE(grid_definition), INTENT(IN) :: outgrid142 REAL(dp), INTENT(IN) :: invar(0:,0:,0:)143 REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:)142 TYPE(grid_definition), INTENT(IN) :: outgrid 143 REAL(dp), INTENT(IN) :: invar(0:,0:,0:) 144 REAL(dp), INTENT(OUT) :: outvar(0:,0:,0:) 144 145 TYPE(nc_var), INTENT(IN), OPTIONAL :: ncvar 145 146 … … 413 414 414 415 END SUBROUTINE rotate_to_cosmo 416 415 417 416 418 … … 427 429 !> ------------- 428 430 !> jj, lat 429 !> ^ j430 !> | \ i431 !> ^ j 432 !> | \ i 431 433 !> jj(i,j,2/3) + ... 2 ---\--------/------ 3 432 434 !> | | ^ \ / | … … 459 461 !> 460 462 !------------------------------------------------------------------------------! 461 SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, cosmo_dxi, & 462 cosmo_dyi, palm_clat, palm_clon, palm_ii, palm_jj) 463 SUBROUTINE find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 464 palm_clat, palm_clon, & 465 palm_ii, palm_jj) 463 466 464 467 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 465 468 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 466 REAL(dp) , INTENT(IN):: cosmo_dxi, cosmo_dyi469 REAL(dp) :: cosmo_dxi, cosmo_dyi 467 470 INTEGER, DIMENSION(0:,0:,1:), INTENT(OUT) :: palm_ii, palm_jj 468 471 … … 472 475 lon0 = cosmo_lon(0) 473 476 lat0 = cosmo_lat(0) 477 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) 478 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) 474 479 475 480 DO j = 0, UBOUND(palm_clon, 2)!palm_grid % ny … … 508 513 TYPE(grid_definition), INTENT(IN) :: palm_intermediate 509 514 510 INTEGER :: i, j, k, nx, ny, nz, nlev, k cur515 INTEGER :: i, j, k, nx, ny, nz, nlev, k_intermediate 511 516 LOGICAL :: point_is_below_grid, point_is_above_grid, & 512 517 point_is_in_current_cell … … 523 528 DO j = 0, ny 524 529 525 k cur= 0530 k_intermediate = 0 526 531 527 532 column_base = palm_intermediate % h(i,j,0) … … 532 537 ! cell, or above column_top. Keep increasing current cell index until 533 538 ! the current cell overlaps with the current_height. 534 DO k = 0, nz539 DO k = 1, nz 535 540 536 541 ! Memorize the top and bottom boundaries of the coarse cell and the 537 542 ! current height within it 538 543 current_height = palm_grid % z(k) + palm_grid % z0 539 h_top = palm_intermediate % h(i,j,k cur+1)540 h_bottom = palm_intermediate % h(i,j,k cur)544 h_top = palm_intermediate % h(i,j,k_intermediate+1) 545 h_bottom = palm_intermediate % h(i,j,k_intermediate) 541 546 542 547 point_is_above_grid = (current_height > column_top) !22000m, very unlikely … … 556 561 palm_grid % w_verti(i,j,k,1:2) = - 2.0_dp 557 562 563 message = "PALM-4U grid extends above COSMO-DE model top." 564 CALL abort('find_vertical_neighbours_and_weights', message) 565 558 566 ELSE IF (point_is_below_grid) THEN 559 567 … … 564 572 ! cycle through intermediate levels until current 565 573 ! intermediate-grid cell overlaps with current_height 566 DO WHILE (.NOT. point_is_in_current_cell .AND. k cur<= nlev-1)567 k cur = kcur+ 1568 569 h_top = palm_intermediate % h(i,j,k cur+1)570 h_bottom = palm_intermediate % h(i,j,k cur)574 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 575 k_intermediate = k_intermediate + 1 576 577 h_top = palm_intermediate % h(i,j,k_intermediate+1) 578 h_bottom = palm_intermediate % h(i,j,k_intermediate) 571 579 point_is_in_current_cell = ( & 572 580 current_height >= h_bottom .AND. & … … 575 583 END DO 576 584 577 ! kcur = 48 indicates the last section (indices 48 and 49), i.e. 578 ! kcur = 49 is not the beginning of a valid cell. 579 IF (kcur > nlev-1) THEN 580 message = "Index " // TRIM(str(kcur)) // " is above intermediate grid range." 585 ! k_intermediate = 48 indicates the last section (indices 48 and 49), i.e. 586 ! k_intermediate = 49 is not the beginning of a valid cell. 587 IF (k_intermediate > nlev-1) THEN 588 message = "Index " // TRIM(str(k_intermediate)) // & 589 " is above intermediate grid range." 581 590 CALL abort('find_vertical_neighbours', message) 582 591 END IF 583 592 584 palm_grid % kk(i,j,k,1) = k cur585 palm_grid % kk(i,j,k,2) = k cur+ 1593 palm_grid % kk(i,j,k,1) = k_intermediate 594 palm_grid % kk(i,j,k,2) = k_intermediate + 1 586 595 587 596 ! copmute vertical weights … … 643 652 ! 644 653 SUBROUTINE compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 645 cosmo_dxi, cosmo_dyi,palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz)654 palm_clat, palm_clon, palm_ii, palm_jj, palm_w_horiz) 646 655 647 656 REAL(dp), DIMENSION(0:), INTENT(IN) :: cosmo_lat, cosmo_lon 648 REAL(dp) , INTENT(IN):: cosmo_dxi, cosmo_dyi657 REAL(dp) :: cosmo_dxi, cosmo_dyi 649 658 REAL(dp), DIMENSION(0:,0:), INTENT(IN) :: palm_clat, palm_clon 650 659 INTEGER, DIMENSION(0:,0:,1:), INTENT(IN) :: palm_ii, palm_jj … … 654 663 REAL(dp) :: wl, wp 655 664 INTEGER :: i, j 665 666 cosmo_dxi = 1.0_dp / (cosmo_lon(1) - cosmo_lon(0)) 667 cosmo_dyi = 1.0_dp / (cosmo_lat(1) - cosmo_lat(0)) 656 668 657 669 DO j = 0, UBOUND(palm_clon, 2) -
palm/trunk/UTIL/inifor/src/types.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Introduced new PALM grid stretching: 24 ! - Converted vertical grid_definition coordinte variables to pointers 25 ! Improved command line interface: 26 ! - Moved INIFOR configuration into a new derived data type 27 ! Removed unnecessary variables 23 28 ! 24 29 ! … … 41 46 42 47 USE defs, & 43 ONLY: dp, PATH, SNAME, LNAME48 ONLY: dp, DATE, PATH, SNAME, LNAME 44 49 USE netcdf, & 45 50 ONLY: NF90_MAX_VAR_DIMS, NF90_MAX_NAME 46 51 47 52 IMPLICIT NONE 53 54 TYPE inifor_config 55 CHARACTER(LEN=DATE) :: start_date !< String of the FORMAT YYYYMMDDHH indicating the start of the intended PALM-4U simulation 56 57 CHARACTER(LEN=PATH) :: input_path !< Path to the input data file directory 58 CHARACTER(LEN=PATH) :: hhl_file !< Path to the file containing the COSMO-DE HHL variable (height of half layers, i.e. vertical cell faces) 59 CHARACTER(LEN=PATH) :: namelist_file !< Path to the PALM-4U namelist file 60 CHARACTER(LEN=PATH) :: output_file !< Path to the INIFOR output file (i.e. PALM-4U dynamic driver') 61 CHARACTER(LEN=PATH) :: soiltyp_file !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types) 62 CHARACTER(LEN=PATH) :: static_driver_file !< Path to the file containing the COSMO-DE SOILTYP variable (map of COSMO-DE soil types) 63 64 CHARACTER(LEN=SNAME) :: flow_prefix !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 65 CHARACTER(LEN=SNAME) :: soil_prefix !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 66 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses 67 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses 68 69 CHARACTER(LEN=SNAME) :: bc_mode 70 CHARACTER(LEN=SNAME) :: ic_mode 71 CHARACTER(LEN=SNAME) :: rotation_method 72 73 REAL(dp) :: p0 74 REAL(dp) :: ug 75 REAL(dp) :: vg 76 REAL(dp) :: z0 !< Elevation of the PALM-4U domain above sea level [m] 77 END TYPE inifor_config 48 78 49 79 TYPE grid_definition … … 55 85 INTEGER, ALLOCATABLE :: jj(:,:,:) !< Given a point (i,j,k) in the PALM-4U grid, jj(i,j,l) gives the y index of the l'th horizontl neighbour on the COSMO-DE grid. 56 86 INTEGER, ALLOCATABLE :: kk(:,:,:,:) !< Given a point (i,j,k) in the PALM-4U grid, kk(i,j,k,l) gives the z index of the l'th vertical neighbour in the intermediate grid. 57 REAL(dp) :: dx !< grid spacing in the first dimension [m]58 REAL(dp) :: dy !< grid spacing in the second dimension [m]59 REAL(dp) :: dz !< grid spacing in the third dimension [m]60 REAL(dp) :: dxi !< inverse grid spacing in the first dimension [m^-1]61 REAL(dp) :: dyi !< inverse grid spacing in the second dimension [m^-1]62 REAL(dp) :: dzi !< inverse grid spacing in the third dimension [m^-1]63 87 REAL(dp) :: lx !< domain length in the first dimension [m] 64 88 REAL(dp) :: ly !< domain length in the second dimension [m] 65 REAL(dp) :: lz !< domain length in the third dimension [m]66 89 REAL(dp) :: x0 !< x coordinate of PALM-4U domain projection centre, i.e. location of zero distortion 67 90 REAL(dp) :: y0 !< y coordinate of PALM-4U domain projection centre, i.e. location of zwro distortion … … 69 92 REAL(dp), ALLOCATABLE :: x(:) !< coordinates of cell centers in x direction [m] 70 93 REAL(dp), ALLOCATABLE :: y(:) !< coordinates of cell centers in y direction [m] 71 REAL(dp), ALLOCATABLE:: z(:) !< coordinates of cell centers in z direction [m]94 REAL(dp), POINTER :: z(:) !< coordinates of cell centers in z direction [m] 72 95 REAL(dp), ALLOCATABLE :: h(:,:,:) !< heights grid point for intermediate grids [m] 73 96 REAL(dp), POINTER :: hhl(:,:,:) !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from … … 76 99 REAL(dp), ALLOCATABLE :: xu(:) !< coordinates of cell faces in x direction [m] 77 100 REAL(dp), ALLOCATABLE :: yv(:) !< coordinates of cell faces in y direction [m] 78 REAL(dp), ALLOCATABLE:: zw(:) !< coordinates of cell faces in z direction [m]101 REAL(dp), POINTER :: zw(:) !< coordinates of cell faces in z direction [m] 79 102 REAL(dp), ALLOCATABLE :: lat(:) !< rotated-pole latitudes of scalars (cell centers) of the COSMO-DE grid [rad] 80 103 REAL(dp), ALLOCATABLE :: lon(:) !< rotated-pole longitudes of scalars (cell centres) of the COSMO-DE grid [rad] … … 89 112 REAL(dp), ALLOCATABLE :: w_horiz(:,:,:) !< weights for bilinear horizontal interpolation 90 113 REAL(dp), ALLOCATABLE :: w_verti(:,:,:,:) !< weights for linear vertical interpolation 91 END TYPE 114 END TYPE grid_definition 92 115 93 116 … … 103 126 INTEGER :: dimvarids_soil(3)!< NetCDF IDs of the grid coordinates for soil points x, y, depth 104 127 REAL(dp), POINTER :: time(:) ! vector of output time steps 105 END TYPE 128 END TYPE nc_file 106 129 107 130 … … 123 146 CHARACTER(LEN=SNAME) :: kind !< Kind of grid 124 147 CHARACTER(LEN=SNAME) :: task !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile' 125 LOGICAL :: to_be_processed = .FALSE. !< I niforflag indicating whether variable shall be processed126 LOGICAL :: is_read = .FALSE. !< I niforflag indicating whether variable has been read127 LOGICAL :: is_upside_down = .FALSE. !< I niforflag indicating whether variable shall be processed148 LOGICAL :: to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed 149 LOGICAL :: is_read = .FALSE. !< INIFOR flag indicating whether variable has been read 150 LOGICAL :: is_upside_down = .FALSE. !< INIFOR flag indicating whether variable shall be processed 128 151 TYPE(grid_definition), POINTER :: grid !< Pointer to the corresponding output grid 129 152 TYPE(grid_definition), POINTER :: intermediate_grid !< Pointer to the corresponding intermediate grid -
palm/trunk/UTIL/inifor/src/util.f90
r2718 r3182 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! Improved real-to-string conversion 23 24 ! 24 25 ! … … 43 44 ONLY : C_CHAR, C_INT, C_PTR, C_SIZE_T 44 45 USE defs, & 45 ONLY : dp, PI, DATE 46 ONLY : dp, PI, DATE, SNAME 46 47 47 48 IMPLICIT NONE … … 279 280 ! Convert a real number to a string in scientific notation 280 281 ! showing four significant digits. 281 CHARACTER(LEN= 11) FUNCTION real_to_str(val, format)282 CHARACTER(LEN=SNAME) FUNCTION real_to_str(val, format) 282 283 283 284 REAL(dp), INTENT(IN) :: val … … 285 286 286 287 IF (PRESENT(format)) THEN 287 WRITE(real_to_str, TRIM(format)) val288 WRITE(real_to_str, format) val 288 289 ELSE 289 290 WRITE(real_to_str, '(E11.4)') val 290 real_to_str = ADJUSTL(real_to_str)291 291 END IF 292 real_to_str = ADJUSTL(real_to_str) 292 293 293 294 END FUNCTION real_to_str
Note: See TracChangeset
for help on using the changeset viewer.
