Changeset 3395
- Timestamp:
- Oct 22, 2018 5:32:49 PM (5 years ago)
- Location:
- palm/trunk/UTIL/inifor
- Files:
-
- 2 deleted
- 12 edited
-
Makefile (modified) (2 diffs)
-
Makefile.gnu (modified) (2 diffs)
-
Makefile.ifort (modified) (2 diffs)
-
README (modified) (12 diffs)
-
src/control.f90 (modified) (3 diffs)
-
src/defs.f90 (modified) (3 diffs)
-
src/get_netcdf_variable.inc (deleted)
-
src/grid.f90 (modified) (94 diffs)
-
src/inifor.f90 (modified) (10 diffs)
-
src/io.f90 (modified) (16 diffs)
-
src/transform.f90 (modified) (11 diffs)
-
src/types.f90 (modified) (7 diffs)
-
src/util.f90 (modified) (4 diffs)
-
tests/test-stretching.f90 (deleted)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/Makefile
r3183 r3395 25 25 # ----------------- 26 26 # $Id$ 27 # Updated build order 28 # 29 # 3183 2018-07-27 14:25:55Z suehring 27 30 # Added __netcdf4 preprocessor flag 28 31 # … … 47 50 TEST_PATH = $(PROJECT_PATH)/tests 48 51 49 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/ util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod52 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod 50 53 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90 51 54 OBJECTS = $(SOURCES:%.f90=%.o) -
palm/trunk/UTIL/inifor/Makefile.gnu
r3183 r3395 25 25 # ----------------- 26 26 # $Id$ 27 # Updated build order 28 # 29 # 3183 2018-07-27 14:25:55Z suehring 27 30 # Added __netcdf4 preprocessor flag 28 31 # Corrected compilation order … … 44 47 PROJECT = inifor 45 48 PROJECT_PATH = . 46 BIN_PATH = $(PROJECT_PATH)/ bin49 BIN_PATH = $(PROJECT_PATH)/../../SCRIPTS 47 50 SRC_PATH = $(PROJECT_PATH)/src 48 51 TEST_PATH = $(PROJECT_PATH)/tests 49 52 50 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/ util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod53 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod 51 54 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90 52 55 OBJECTS = $(SOURCES:%.f90=%.o) -
palm/trunk/UTIL/inifor/Makefile.ifort
r3183 r3395 25 25 # ----------------- 26 26 # $Id$ 27 # Updated build order 28 # 29 # 3183 2018-07-27 14:25:55Z suehring 27 30 # Added __netcdf4 preprocessor flag 28 31 # Corrected compilation order … … 44 47 PROJECT = inifor 45 48 PROJECT_PATH = . 46 BIN_PATH = $(PROJECT_PATH)/ bin49 BIN_PATH = $(PROJECT_PATH)/../../SCRIPTS 47 50 SRC_PATH = $(PROJECT_PATH)/src 48 51 TEST_PATH = $(PROJECT_PATH)/tests 49 52 50 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/ util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod53 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod 51 54 SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90 52 55 OBJECTS = $(SOURCES:%.f90=%.o) -
palm/trunk/UTIL/inifor/README
r3182 r3395 1 # INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1. 3.0)1 # INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.0) 2 2 3 3 INIFOR provides the meteorological fields required to initialize and drive the … … 14 14 15 15 16 ## MINIMAL USAGE EXAMPLE 17 18 1. Customize `./namelist` (number of grid points and spacings, end_time) 19 2. Run `current_version/trunk/SCRIPTS/inifor --path <scenario path> --date <YYYYMMDD>` 20 21 16 22 ## USAGE 17 23 18 1. Customize `./namelist` (number or grid points and spacings, end_time) 19 2. Run `current_version/trunk/SCRIPTS/inifor -path <scenario path> -date <YYYYMMDD>` 24 After compilation, the `inifor` binary resides in the `$PALM_BIN` path, i.e. in 25 `<PALM path>/current_version/trunk/SCRIPTS/`. 26 27 In order to run, INIFOR requires three kinds of inputs: 28 29 1. hourly COSMO model output, 30 2. a steering namelist file, and 31 3. command-line options. 32 33 In addition, a static driver file may be supplied (see `--static-driver` option 34 below) in order to pass the coordinates of the PALM origin to INIFOR. If no 35 static driver is passed to INIFOR, origin coordinates are read from the namelist 36 file. 37 38 A typical `inifor` call looks like this: 39 40 inifor --path /data/evaluation/20170729 --date 2017073006 \ 41 --init-mode profile -n namelist -o dynamic_driver.nc \ 42 --elevation 110 --input-prefix lff0 43 44 ### INPUT DATA: COSMO MODEL OUTPUT 45 46 INIFOR processes COSMO model output which it requires to be stored in a set of 47 netCDF files located in a user-specified path (see `--path` option). These are: 48 49 - hhl.nc: This file provides the COSMO numerical grid. (*hhl* abbreviates 50 *height of half layers*, i.e. the heights of the vertical cell boundaries.) 51 - soil.nc: This file provides the COSMO soil map and is used to destinguish land 52 from water cells. 53 - `<prefix>YYYYMMDDHH-<suffix>.nc`: Each of these files contains COSMO model 54 of one time step at the given time in UTC (Y..year, M..month, D..day, H..hour). 55 - The `<prefix>` distinguishes different DWD products, for instance COSMO 56 analyses (`laf`) or forecasts (`lff`). 57 - The `<suffix>` distinguishes four kinds of COSMO model output data, namely 58 - `flow` (atmospheric fields) 59 - `rad` (radiation) 60 - `soil` (soil moisture and temperature) 61 - `soilmoisture` (precipitation and evaporation) 62 - For example, laf2016010100-flow.nc contains the atmospheric fields of the 63 COSMO analysis of Januray 1st, 2016 for 0:00 UTC. 20 64 21 65 22 ## AVAILABLE NAMELIST PARAMETERS66 ### AVAILABLE NAMELIST PARAMETERS 23 67 24 68 INIFOR mirrors a subset of the PALM-4U's Fortran namelists `inipar` and `d3par` … … 26 70 27 71 28 ### inipar72 #### inipar 29 73 30 74 nx, ny, nz - number of PALM-4U grid points in x, y, and z direction … … 34 78 dz_max - maximum vertical grid spacing [m] 35 79 dz_stretch_level_start(9) - array of height levels above which the grid is 36 to be stretched vertically [m]80 to be stretched vertically [m] 37 81 dz_stretch_level_end(9) - array of height levels until which the grid is to 38 be stretched vertically [m]82 be stretched vertically [m] 39 83 longitude, latitude - geographical coordinates of the PALM-4U origin [deg] 40 84 41 85 42 ### d3par86 #### d3par 43 87 44 88 end_time - PALM-4U simulation time. INIFOR will produce hourly forcing data … … 47 91 48 92 49 ### EXAMPLE NAMELIST FILE93 #### EXAMPLE NAMELIST FILE 50 94 51 95 &inipar nx = 4679, ny = 3939, nz = 360 … … 59 103 60 104 61 ## AVAILABLE COMMAND-LINE PARAMETERS 105 ### AVAILABLE COMMAND-LINE PARAMETERS 106 --averaging-mode <mode>: 107 Selects how averaged quantities (large-scale forcing terms) are computed. 108 INIFOR supports averaging along input model levels ('level') and along 109 constant heights ('height'). Default: level 110 111 -a, --averaging-angle <angle>: 112 Width of the averaging box in longitudal and latitudal direction in the 113 source coordinate system (COSMO rotated-pole) [deg]. Default: 2.0 62 114 63 115 -d, --date, -date <date>: … … 65 117 hours (HH) are given, INIFOR assumes that the simulation starts at O UTC 66 118 on that day. Default: 20130721 67 119 68 120 -i, --init-mode, -mode <mode>: 69 121 Set the PALM-4U initialization mode. INIFOR can provide initial conditions … … 72 124 73 125 -l, --hhl-file, -hhl <netCDF file>: 74 Location of the netCDF file containing the vertical COSMO -DEgrid levels,126 Location of the netCDF file containing the vertical COSMO grid levels, 75 127 specifically the heights of half levels (hhl, i.e. vertical cell faces). 76 128 Default: <scenario path>/hhl.nc … … 83 135 -o, --output <output file>: 84 136 Name of the INIFOR output file, i.e. the PALM-4U dynamic driver. 85 Default: ./ palm-4u-input.nc137 Default: ./dynamic_driver.nc 86 138 87 139 -p, --path, -path <scenario path>: … … 89 141 90 142 -r, --surface-pressure, -p0 <pressure>: 91 Surface pressure at z=0 in the PALM-4U domain [Pa].92 Default: 1e5143 Manually set the pressure at z=0 in the PALM-4U domain [Pa]. If not 144 given, the surface pressure is computed from the COSMO pressure field. 93 145 94 146 -s, --soil-file, -soil <netCDF file>: 95 Location of the netCDF file containing the COSMO -DEsoil type map.147 Location of the netCDF file containing the COSMO soil type map. 96 148 Default: <scenario path>/soil.nc 97 149 98 150 -t, --static-driver, -static <netCDF file>: 99 Location of the netCDF file containing the static driver f ile for the case151 Location of the netCDF file containing the static driver for the case 100 152 to be simulated with PALM-4U. Optional parameter. Default: None 101 153 102 154 -u, --geostrophic-u, -ug <velocity>: 103 Specifies the geostrophic wind in x direction [m/s]. Default: 0 155 Manually specify the geostrophic wind in x direction [m/s]. If not 156 given, the geostrophic wind is computed from the COSMO pressure field. 104 157 105 158 -v, --geostrophic-v, -vg <velocity>: 106 Specifies the geostrophic wind in y direction [m/s]. Default: 0 107 108 --version: 109 Output version number and exit. 159 Manually specify the geostrophic wind in y direction [m/s]. If not 160 given, the geostrophic wind is computed from the COSMO pressure field. 110 161 111 162 -z, --elevation, -z0 <height>: Specifies the elevation of the PALM-4U domain … … 113 164 114 165 115 ## ADDITIONAL COMMAND-LINE OPTIONS 166 #### ADDITIONAL COMMAND-LINE OPTIONS 167 168 --input-prefix <prefix>: 169 Set the file prefixes for all input files. Individual prefixes can be 170 overwritten with the options below. Default: laf 116 171 117 172 --flow-prefix <prefix>: … … 126 181 --soilmoisture-prefix <prefix>: 127 182 Set the file prefix of soil moisture input files. Default: laf 183 184 --debug: 185 Enable debugging messages. 186 187 --version: 188 Output version number and exit. 189 -
palm/trunk/UTIL/inifor/src/control.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Suppress debugging messages unless --debug option is given 29 ! 30 ! 31 ! 3183 2018-07-27 14:25:55Z suehring 28 32 ! Added version and copyright output 29 33 ! … … 57 61 CONTAINS 58 62 59 SUBROUTINE report(routine, message) 60 61 CHARACTER(LEN=*), INTENT(IN) :: routine 62 CHARACTER(LEN=*), INTENT(IN) :: message 63 INTEGER :: u 64 LOGICAL, SAVE :: is_first_run = .TRUE. 65 66 PRINT *, "inifor: " // TRIM(message) // " [ " // TRIM(routine) // " ]" 63 SUBROUTINE report(routine, message, debug) 64 65 CHARACTER(LEN=*), INTENT(IN) :: routine 66 CHARACTER(LEN=*), INTENT(IN) :: message 67 LOGICAL, OPTIONAL, INTENT(IN) :: debug 68 INTEGER :: u 69 LOGICAL, SAVE :: is_first_run = .TRUE. 70 LOGICAL :: suppress_message 71 67 72 68 73 IF (is_first_run) THEN … … 73 78 END IF 74 79 75 WRITE(u, *) TRIM(message) // " [ " // TRIM(routine) // " ]" 80 81 suppress_message = .FALSE. 82 IF (PRESENT(debug)) THEN 83 IF (.NOT. debug) suppress_message = .TRUE. 84 END IF 85 86 IF (.NOT. suppress_message) THEN 87 PRINT *, "inifor: " // TRIM(message) // " [ " // TRIM(routine) // " ]" 88 WRITE(u, *) TRIM(message) // " [ " // TRIM(routine) // " ]" 89 END IF 76 90 77 91 CLOSE(u) -
palm/trunk/UTIL/inifor/src/defs.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! New parameters for computation of geostrophic winds 29 ! Bumped INIFOR version number 30 ! 31 ! 32 ! 3183 2018-07-27 14:25:55Z suehring 28 33 ! Updated defaults for soil extrapolation steps and nudging time-scale 29 34 ! Improved handling of the start date string … … 73 78 REAL(dp), PARAMETER :: G = 9.80665_dp !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s] 74 79 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] 80 REAL(dp), PARAMETER :: HECTO = 100_dp !< unit conversion factor from hPa to Pa 75 81 76 82 ! PALM-4U parameters 83 REAL(dp), PARAMETER :: OMEGA = 7.29e-5_dp !< angular velocity of Earth's rotation [s^-1] 77 84 REAL(dp), PARAMETER :: P_REF = 1e5_dp !< Reference pressure for potential temperature [Pa] 78 85 REAL(dp), PARAMETER :: RD_PALM = 287.0_dp !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K] … … 83 90 INTEGER, PARAMETER :: FORCING_STEP = 1 !< Number of hours between forcing time steps [h] 84 91 REAL(dp), PARAMETER :: NUDGING_TAU = 21600.0_dp !< Nudging relaxation time scale [s] 85 CHARACTER(LEN=*), PARAMETER :: VERSION = '1. 3.0' !< INIFOR version number92 CHARACTER(LEN=*), PARAMETER :: VERSION = '1.4.0' !< INIFOR version number 86 93 CHARACTER(LEN=*), PARAMETER :: COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // & 87 94 NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice -
palm/trunk/UTIL/inifor/src/grid.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Added computation of geostrophic winds form COSMO pressure fields 29 ! Introduced averaging grids and internal 'output' variables for computation of 30 ! geostrophic and large-scale forcings 31 ! Merged qv, uv and vg output variables and 'temperature' IO group into new 32 ! 'thermodynamics' IP group 33 ! Cleaned up terminal output, show some messages only with --debug option 34 ! Improved variable naming and readability 35 ! 36 ! 37 ! 3183 2018-07-27 14:25:55Z suehring 28 38 ! Introduced new PALM grid stretching 29 39 ! Updated variable names and metadata for PIDS v1.9 compatibility … … 61 71 ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp, & 62 72 SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, & 63 BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L 73 BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM, & 74 RHO_L, OMEGA, HECTO 64 75 USE io, & 65 76 ONLY: get_netcdf_variable, get_netcdf_attribute, & … … 68 79 ONLY: NF90_MAX_NAME, NF90_MAX_VAR_DIMS 69 80 USE transform, & 70 ONLY: rotate_to_cosmo, find_horizontal_neighbours,&81 ONLY: average_2d, rotate_to_cosmo, find_horizontal_neighbours,& 71 82 compute_horizontal_interp_weights, & 72 find_vertical_neighbours_and_weights, interpolate_2d, & 73 gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d, & 74 project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot 83 find_vertical_neighbours_and_weights_interp, & 84 find_vertical_neighbours_and_weights_average, interpolate_2d, & 85 gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, & 86 project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot, & 87 phirot2phi, rlarot2rla 75 88 USE types 76 89 USE util … … 80 93 SAVE 81 94 95 REAL(dp) :: averaging_angle = 0.0_dp !< latitudal and longitudal width of averaging regions [rad] 96 REAL(dp) :: averaging_width_ns = 0.0_dp !< longitudal width of averaging regions [m] 97 REAL(dp) :: averaging_width_ew = 0.0_dp !< latitudal width of averaging regions [m] 82 98 REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad] 83 99 REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad] … … 87 103 REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] 88 104 REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] 105 REAL(dp) :: lam_centre = 0.0_dp !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 106 REAL(dp) :: phi_centre = 0.0_dp !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 107 REAL(dp) :: lam_east = 0.0_dp !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 108 REAL(dp) :: lam_west = 0.0_dp !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 109 REAL(dp) :: phi_north = 0.0_dp !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 110 REAL(dp) :: phi_south = 0.0_dp !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 89 111 REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug 90 112 REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m] … … 100 122 REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1] 101 123 REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1] 124 REAL(dp) :: f3 = 0.0_dp !< Coriolis parameter 102 125 REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m] 103 126 REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m] 104 REAL(dp) :: ug = 0.0_dp !< geostrophic wind in x direction [m/s]105 REAL(dp) :: vg = 0.0_dp !< geostrophic wind in y direction [m/s]106 127 REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa] 107 128 REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m] … … 140 161 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 141 162 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 142 INTEGER :: i !< indexing variable143 INTEGER :: average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging144 INTEGER :: k !< indexing variable145 163 INTEGER :: nt !< number of output time steps 146 164 INTEGER :: nx !< number of PALM-4U grid points in x direction … … 164 182 LOGICAL :: profile_grids_required 165 183 166 INTEGER :: n_invar = 0 !< number of variables in the input variable table167 INTEGER :: n_outvar = 0 !< number of variables in the output variable table168 184 TYPE(nc_var), ALLOCATABLE, TARGET :: input_var_table(:) !< table of input variables 169 185 TYPE(nc_var), ALLOCATABLE, TARGET :: output_var_table(:) !< table of input variables … … 219 235 TYPE(grid_definition), TARGET :: w_south_intermediate !< 220 236 TYPE(grid_definition), TARGET :: w_top_intermediate !< 221 TYPE(grid_definition), TARGET :: scalar_profile_grid !< 222 TYPE(grid_definition), TARGET :: scalar_profile_intermediate !< 223 TYPE(grid_definition), TARGET :: w_profile_grid !< 224 TYPE(grid_definition), TARGET :: w_profile_intermediate !< 237 238 TYPE(grid_definition), TARGET :: north_averaged_scalar_profile 239 TYPE(grid_definition), TARGET :: south_averaged_scalar_profile 240 TYPE(grid_definition), TARGET :: west_averaged_scalar_profile 241 TYPE(grid_definition), TARGET :: east_averaged_scalar_profile 242 TYPE(grid_definition), TARGET :: averaged_scalar_profile 243 TYPE(grid_definition), TARGET :: averaged_w_profile 225 244 226 245 TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable … … 237 256 CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files 238 257 CHARACTER(LEN=SNAME) :: input_prefix !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses 258 CHARACTER(LEN=SNAME) :: flow_prefix !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 259 CHARACTER(LEN=SNAME) :: soil_prefix !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 260 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses 261 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses 239 262 CHARACTER(LEN=SNAME) :: flow_suffix !< Suffix of flow input files, e.g. 'flow' 240 263 CHARACTER(LEN=SNAME) :: soil_suffix !< Suffix of soil input files, e.g. 'soil' … … 249 272 250 273 SUBROUTINE setup_parameters() 274 INTEGER :: k 251 275 252 276 ! … … 282 306 283 307 ! Default atmospheric parameters 284 ug = 0.0_dp285 vg = 0.0_dp308 cfg % ug = 0.0_dp 309 cfg % vg = 0.0_dp 286 310 cfg % p0 = P_SL 287 311 … … 296 320 input_prefix = 'laf' ! 'laf' for COSMO-DE analyses 297 321 cfg % flow_prefix = input_prefix 322 cfg % input_prefix = input_prefix 298 323 cfg % soil_prefix = input_prefix 299 324 cfg % radiation_prefix = input_prefix … … 304 329 radiation_suffix = '-rad' 305 330 soilmoisture_suffix = '-soilmoisture' 331 332 cfg % debug = .FALSE. 333 cfg % averaging_angle = 2.0_dp 306 334 ! 307 335 !------------------------------------------------------------------------------ … … 318 346 cfg % ic_mode = 'volume' 319 347 cfg % bc_mode = 'real' 320 321 ! Update default file names and domain centre 348 cfg % averaging_mode = 'level' 349 350 ! Overwrite defaults with user configuration 322 351 CALL parse_command_line_arguments( cfg ) 352 353 flow_prefix = TRIM(cfg % input_prefix) 354 radiation_prefix = TRIM(cfg % input_prefix) 355 soil_prefix = TRIM(cfg % input_prefix) 356 soilmoisture_prefix = TRIM(cfg % input_prefix) 357 IF (cfg % flow_prefix_is_set) flow_prefix = TRIM(cfg % flow_prefix) 358 IF (cfg % radiation_prefix_is_set) radiation_prefix = TRIM(cfg % radiation_prefix) 359 IF (cfg % soil_prefix_is_set) soil_prefix = TRIM(cfg % soil_prefix) 360 IF (cfg % soilmoisture_prefix_is_set) soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix) 323 361 324 362 output_file % name = cfg % output_file … … 336 374 END IF 337 375 376 377 ! Set default file paths, if not specified by user. 338 378 CALL normalize_path(cfg % input_path) 339 379 IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' … … 361 401 362 402 ! Generate input file lists 363 CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour, & 364 cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files) 365 CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour, & 366 cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files) 367 CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, & 368 cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files) 369 CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 370 cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) 403 CALL get_input_file_list( & 404 cfg % start_date, start_hour_flow, end_hour, step_hour, & 405 cfg % input_path, flow_prefix, flow_suffix, flow_files) 406 CALL get_input_file_list( & 407 cfg % start_date, start_hour_soil, end_hour, step_hour, & 408 cfg % input_path, soil_prefix, soil_suffix, soil_files) 409 CALL get_input_file_list( & 410 cfg % start_date, start_hour_radiation, end_hour, step_hour, & 411 cfg % input_path, radiation_prefix, radiation_suffix, radiation_files) 412 CALL get_input_file_list( & 413 cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 414 cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) 371 415 372 416 ! … … 499 543 500 544 CALL run_control('time', 'comp') 545 546 !------------------------------------------------------------------------------ 547 ! Section 4.3: INIFOR averaging domains 548 !------------------------------------------------------------------------------ 549 550 ! Compute coordiantes of the PALM centre in the source (COSMO) system 551 phi_centre = phirot2phi( & 552 phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 553 rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 554 polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, & 555 polgam = gam * TO_DEGREES & 556 ) * TO_RADIANS 557 558 lam_centre = rlarot2rla( & 559 phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 560 rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 561 polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, & 562 polgam = gam * TO_DEGREES & 563 ) * TO_RADIANS 564 565 message = "PALM-4U centre:" // NEW_LINE('') // & 566 " lon (lambda) = " // & 567 TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & 568 " lat (phi ) = " // & 569 TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" 570 CALL report( 'setup_parameters', message ) 571 572 573 ! Compute boundaries of the central averaging box 574 averaging_angle = cfg % averaging_angle * TO_RADIANS 575 lam_east = lam_centre + 0.5_dp * averaging_angle 576 lam_west = lam_centre - 0.5_dp * averaging_angle 577 phi_north = phi_centre + 0.5_dp * averaging_angle 578 phi_south = phi_centre - 0.5_dp * averaging_angle 579 averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS 580 averaging_width_ns = averaging_angle * EARTH_RADIUS 581 582 ! Coriolis parameter 583 f3 = 2.0_dp * OMEGA * SIN( & 584 TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,& 585 phi_n * TO_DEGREES, lambda_n * TO_DEGREES, & 586 gam * TO_DEGREES ) & 587 ) 501 588 502 589 END SUBROUTINE setup_parameters … … 603 690 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 604 691 x0=x0, y0=y0, z0 = z0, & 605 nx = nx, ny = ny, nz = nlev - 2)692 nx = nx, ny = ny, nz = nlev - 1) 606 693 607 694 IF (boundary_variables_required) THEN … … 824 911 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 825 912 x0=x0, y0=y0, z0 = z0, & 826 nx = 0, ny = ny, nz = nlev - 2)913 nx = 0, ny = ny, nz = nlev - 1) 827 914 828 915 CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & … … 830 917 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 831 918 x0=x0, y0=y0, z0 = z0, & 832 nx = 0, ny = ny, nz = nlev - 2)919 nx = 0, ny = ny, nz = nlev - 1) 833 920 834 921 CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & … … 836 923 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 837 924 x0=x0, y0=y0, z0 = z0, & 838 nx = nx, ny = 0, nz = nlev - 2)925 nx = nx, ny = 0, nz = nlev - 1) 839 926 840 927 CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & … … 842 929 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 843 930 x0=x0, y0=y0, z0 = z0, & 844 nx = nx, ny = 0, nz = nlev - 2)931 nx = nx, ny = 0, nz = nlev - 1) 845 932 846 933 CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & … … 848 935 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 849 936 x0=x0, y0=y0, z0 = z0, & 850 nx = nx, ny = ny, nz = nlev - 2)937 nx = nx, ny = ny, nz = nlev - 1) 851 938 END IF 852 939 … … 856 943 !------------------------------------------------------------------------------ 857 944 858 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & 945 CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid, & 946 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 947 lonmin = lam_west, lonmax = lam_east, & 948 latmin = phi_south, latmax = phi_north, & 949 kind='scalar') 950 951 CALL init_averaging_grid(averaged_w_profile, cosmo_grid, & 952 x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, & 953 lonmin = lam_west, lonmax = lam_east, & 954 latmin = phi_south, latmax = phi_north, & 955 kind='w') 956 957 CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid, & 958 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 959 lonmin = lam_west, lonmax = lam_east, & 960 latmin = phi_centre - averaging_angle, & 961 latmax = phi_centre, & 962 kind='scalar') 963 964 CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid, & 965 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 966 lonmin = lam_west, lonmax = lam_east, & 967 latmin = phi_centre, & 968 latmax = phi_centre + averaging_angle, & 969 kind='scalar') 970 971 CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid, & 972 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 973 lonmin = lam_centre - averaging_angle, & 974 lonmax = lam_centre, & 975 latmin = phi_south, latmax = phi_north, & 976 kind='scalar') 977 978 CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid, & 979 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, & 980 lonmin = lam_centre, & 981 lonmax = lam_centre + averaging_angle, & 982 latmin = phi_south, latmax = phi_north, & 983 kind='scalar') 984 985 profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & 859 986 TRIM(cfg % bc_mode) == 'ideal' ) 860 987 861 IF (profile_grids_required) THEN 862 CALL init_grid_definition('boundary', grid=scalar_profile_grid, & 863 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 864 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 865 x0=x0, y0=y0, z0 = z0, & 866 nx = 0, ny = 0, nz = nz, z=z) 867 868 CALL init_grid_definition('boundary', grid=w_profile_grid, & 869 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 870 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 871 x0=x0, y0=y0, z0 = z0, & 872 nx = 0, ny = 0, nz = nz - 1, z=zw) 873 874 CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& 875 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 876 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 877 x0=x0, y0=y0, z0 = z0, & 878 nx = 0, ny = 0, nz = nlev - 2, z=z) 879 880 CALL init_grid_definition('boundary', grid=w_profile_intermediate, & 881 xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 882 ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 883 x0=x0, y0=y0, z0 = z0, & 884 nx = 0, ny = 0, nz = nlev - 2, z=zw) 885 END IF 988 !IF (profile_grids_required) THEN 989 ! ! Here, I use the 'boundary' kind, so the vertical coordinate uses the 990 ! ! PALM z coordinate. 991 ! CALL init_grid_definition('boundary', grid=scalar_profile_grid, & 992 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 993 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 994 ! x0=x0, y0=y0, z0 = z0, & 995 ! nx = 0, ny = 0, nz = nz, z=z) 996 997 ! CALL init_grid_definition('boundary', grid=w_profile_grid, & 998 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 999 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1000 ! x0=x0, y0=y0, z0 = z0, & 1001 ! nx = 0, ny = 0, nz = nz - 1, z=zw) 1002 1003 ! CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& 1004 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 1005 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1006 ! x0=x0, y0=y0, z0 = z0, & 1007 ! nx = 0, ny = 0, nz = nlev - 2, z=z) 1008 1009 ! CALL init_grid_definition('boundary', grid=w_profile_intermediate, & 1010 ! xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & 1011 ! ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & 1012 ! x0=x0, y0=y0, z0 = z0, & 1013 ! nx = 0, ny = 0, nz = nlev - 2, z=zw) 1014 !END IF 886 1015 887 1016 ! … … 930 1059 931 1060 IF (TRIM(cfg % ic_mode) == 'profile') THEN 932 CALL setup_averaging(cosmo_grid, palm_intermediate, & 933 average_imin, average_imax, average_jmin, average_jmax) 1061 !TODO: remove this conditional if not needed. 934 1062 END IF 935 1063 … … 938 1066 939 1067 940 SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)1068 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode) 941 1069 942 1070 TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid 943 TYPE(grid_definition), INTENT(INOUT), TARGET :: palm_grid, palm_intermediate1071 TYPE(grid_definition), INTENT(INOUT), TARGET :: grid, intermediate_grid 944 1072 CHARACTER, INTENT(IN) :: kind 945 1073 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 946 1074 947 TYPE(grid_definition), POINTER :: grid 948 REAL(dp), DIMENSION(:), POINTER :: lat, lon 949 REAL(dp), DIMENSION(:,:,:), POINTER :: h 950 951 LOGICAL :: setup_vertical 952 953 setup_vertical = .TRUE. 954 IF (PRESENT(ic_mode)) THEN 955 IF (TRIM(ic_mode) == 'profile') setup_vertical = .FALSE. 956 END IF 1075 REAL(dp), DIMENSION(:), POINTER :: cosmo_lat, cosmo_lon 1076 REAL(dp), DIMENSION(:,:,:), POINTER :: cosmo_h 1077 1078 LOGICAL :: setup_volumetric 957 1079 958 1080 !------------------------------------------------------------------------------ … … 964 1086 CASE('s') ! scalars 965 1087 966 lat => cosmo_grid % lat967 lon => cosmo_grid % lon968 h => cosmo_grid % hfl1088 cosmo_lat => cosmo_grid % lat 1089 cosmo_lon => cosmo_grid % lon 1090 cosmo_h => cosmo_grid % hfl 969 1091 970 1092 CASE('w') ! vertical velocity 971 1093 972 lat => cosmo_grid % lat973 lon => cosmo_grid % lon974 h => cosmo_grid % hhl1094 cosmo_lat => cosmo_grid % lat 1095 cosmo_lon => cosmo_grid % lon 1096 cosmo_h => cosmo_grid % hhl 975 1097 976 1098 CASE('u') ! x velocity 977 1099 978 lat => cosmo_grid % lat979 lon => cosmo_grid % lonu980 h => cosmo_grid % hfl1100 cosmo_lat => cosmo_grid % lat 1101 cosmo_lon => cosmo_grid % lonu 1102 cosmo_h => cosmo_grid % hfl 981 1103 982 1104 CASE('v') ! y velocity 983 1105 984 lat => cosmo_grid % latv985 lon => cosmo_grid % lon986 h => cosmo_grid % hfl1106 cosmo_lat => cosmo_grid % latv 1107 cosmo_lon => cosmo_grid % lon 1108 cosmo_h => cosmo_grid % hfl 987 1109 988 1110 CASE DEFAULT … … 993 1115 END SELECT 994 1116 995 grid => palm_intermediate 996 997 CALL find_horizontal_neighbours(lat, lon, & 998 grid % clat, grid % clon, grid % ii, grid % jj) 999 1000 CALL compute_horizontal_interp_weights(lat, lon, & 1001 grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz) 1117 CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 1118 intermediate_grid % clat, intermediate_grid % clon, & 1119 intermediate_grid % ii, intermediate_grid % jj) 1120 1121 CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 1122 intermediate_grid % clat, intermediate_grid % clon, & 1123 intermediate_grid % ii, intermediate_grid % jj, & 1124 intermediate_grid % w_horiz) 1002 1125 1003 1126 !------------------------------------------------------------------------------ … … 1005 1128 !------------------------------------------------------------------------------ 1006 1129 1007 IF (setup_vertical) THEN 1008 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 1009 grid % h(:,:,:) = - EARTH_RADIUS 1130 ! If profile initialization is chosen, we--somewhat counterintuitively-- 1131 ! don't need to compute vertical interpolation weights. At least, we 1132 ! don't need them on the intermediate grid, which fills the entire PALM 1133 ! domain volume. Instead we need vertical weights for the intermediate 1134 ! profile grids, which get computed in setup_averaging(). 1135 setup_volumetric = .TRUE. 1136 IF (PRESENT(ic_mode)) THEN 1137 IF (TRIM(ic_mode) == 'profile') setup_volumetric = .FALSE. 1138 END IF 1139 1140 IF (setup_volumetric) THEN 1141 ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx, & 1142 0:intermediate_grid % ny, & 1143 0:intermediate_grid % nz) ) 1144 intermediate_grid % h(:,:,:) = - EARTH_RADIUS 1010 1145 1011 1146 ! For w points, use hhl, for scalars use hfl 1012 1147 ! compute the full heights for the intermediate grids 1013 CALL interpolate_2d(cosmo_ grid % hfl, grid % h,grid)1014 CALL find_vertical_neighbours_and_weights (palm_grid,grid)1148 CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid) 1149 CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid) 1015 1150 END IF 1016 1151 1017 1152 END SUBROUTINE setup_interpolation 1018 1019 1020 SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)1021 1022 TYPE(grid_definition), INTENT(IN) :: cosmo_grid, palm_intermediate1023 INTEGER, INTENT(INOUT) :: imin, imax, jmin, jmax1024 1025 TYPE(grid_definition), POINTER :: grid1026 REAL(dp) :: lonmin_pos,lonmax_pos, latmin_pos, latmax_pos1027 REAL(dp) :: cosmo_dxi, cosmo_dyi1028 1029 cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))1030 cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))1031 1032 ! find horizontal index ranges for profile averaging1033 lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi1034 lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi1035 latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi1036 latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi1037 1038 imin = FLOOR(lonmin_pos)1039 imax = CEILING(lonmax_pos)1040 jmin = FLOOR(latmin_pos)1041 jmax = CEILING(latmax_pos)1042 1043 ! average heights for intermediate scalar and w profile grids1044 grid => scalar_profile_intermediate1045 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )1046 grid % h(:,:,:) = - EARTH_RADIUS1047 CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax)1048 CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid)1049 1050 grid => w_profile_intermediate1051 ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )1052 grid % h(:,:,:) = - EARTH_RADIUS1053 CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax)1054 CALL find_vertical_neighbours_and_weights(w_profile_grid, grid)1055 1056 END SUBROUTINE setup_averaging1057 1058 1153 1059 1154 !------------------------------------------------------------------------------! … … 1291 1386 1292 1387 END SUBROUTINE init_grid_definition 1388 1389 1390 !------------------------------------------------------------------------------! 1391 ! Description: 1392 ! ------------ 1393 !> Initializes averagin_grid-type variables. 1394 !> 1395 !> Input parameters: 1396 !> ----------------- 1397 !> 1398 !> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid 1399 !> providing longitudes and latitudes of the parent grid 1400 !> 1401 !> x, y : location of the profile to be averaged in the PALM system [m] 1402 !> 1403 !> z : vertical grid coordinates of the profile in the PALM system [m] 1404 !> 1405 !> <lat/lon><min/max>: boundaries of the averaging region in the parent system, 1406 !> i.e. the COSMO-DE rotated-pole system. [rad] 1407 !> 1408 !> kind : kind of quantity to be averaged using this averaging grid. 1409 !> Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the 1410 !> parent/COSMO columns for the region in get_cosmo_averaging_region() is 1411 !> independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity 1412 !> components, which are computed in the preprocessing step. 1413 !> 1414 !> Output parameters: 1415 !> ------------------ 1416 !> avg_grid : averaging_grid variable to be initialized. 1417 !------------------------------------------------------------------------------! 1418 SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0, & 1419 lonmin, lonmax, latmin, latmax, kind) 1420 1421 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1422 TYPE(grid_definition), INTENT(IN) :: cosmo_grid 1423 REAL(dp), INTENT(IN) :: x, y, z0 1424 REAL(dp), INTENT(IN), TARGET :: z(:) 1425 REAL(dp), INTENT(IN) :: lonmin, lonmax, latmin, latmax 1426 1427 CHARACTER(LEN=*), INTENT(IN) :: kind 1428 1429 LOGICAL :: level_based_averaging 1430 1431 INTEGER :: i, j 1432 1433 ALLOCATE( avg_grid % x(1) ) 1434 ALLOCATE( avg_grid % y(1) ) 1435 avg_grid % x(1) = x 1436 avg_grid % y(1) = y 1437 avg_grid % z => z 1438 avg_grid % z0 = z0 1439 1440 avg_grid % nz = SIZE(z, 1) 1441 1442 ALLOCATE( avg_grid % lon(2) ) 1443 ALLOCATE( avg_grid % lat(2) ) 1444 avg_grid % lon(1:2) = (/lonmin, lonmax/) 1445 avg_grid % lat(1:2) = (/latmin, latmax/) 1446 1447 avg_grid % kind = TRIM(kind) 1448 1449 ! Find and store COSMO columns that fall into the coordinate range 1450 ! given by avg_grid % clon, % clat 1451 CALL get_cosmo_averaging_region(avg_grid, cosmo_grid) 1452 1453 ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) ) 1454 ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) ) 1455 ! Compute average COSMO levels in the averaging region 1456 SELECT CASE(avg_grid % kind) 1457 1458 CASE('scalar', 'u', 'v') 1459 avg_grid % cosmo_h => cosmo_grid % hfl 1460 1461 CASE('w') 1462 avg_grid % cosmo_h => cosmo_grid % hhl 1463 1464 CASE DEFAULT 1465 message = "Averaging grid kind '" // TRIM(avg_grid % kind) // & 1466 "' is not supported. Use 'scalar', 'u', or 'v'." 1467 CALL abort('get_cosmo_averaging_region', message) 1468 1469 END SELECT 1470 1471 ! For level-besed averaging, compute average heights 1472 level_based_averaging = .TRUE. 1473 IF (level_based_averaging) THEN 1474 ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) ) 1475 1476 CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:), & 1477 avg_grid % iii, avg_grid % jjj) 1478 END IF 1479 1480 ! Copy averaged grid to all COSMO columnts, leads to computing the same 1481 ! vertical interpolation weights for all columns and to COSMO grid level 1482 ! based averaging onto the averaged COSMO heights. 1483 IF ( TRIM(cfg % averaging_mode) == 'level' ) THEN 1484 FORALL( & 1485 i=1:SIZE(avg_grid % cosmo_h, 1), & 1486 j=1:SIZE(avg_grid % cosmo_h, 2) & 1487 ) avg_grid % cosmo_h(i,j,:) = avg_grid % h(1,1,:) 1488 END IF 1489 1490 ! Compute vertical weights and neighbours 1491 CALL find_vertical_neighbours_and_weights_average(avg_grid) 1492 1493 END SUBROUTINE init_averaging_grid 1494 1495 1496 SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid) 1497 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1498 TYPE(grid_definition), TARGET, INTENT(IN) :: cosmo_grid 1499 1500 REAL(dp), DIMENSION(:), POINTER :: cosmo_lon, cosmo_lat 1501 REAL(dp) :: dlon, dlat 1502 1503 INTEGER :: i, j, imin, imax, jmin, jmax, l, nx, ny 1504 1505 1506 SELECT CASE( TRIM(avg_grid % kind) ) 1507 1508 CASE('scalar', 'w') 1509 cosmo_lon => cosmo_grid % lon 1510 cosmo_lat => cosmo_grid % lat 1511 1512 CASE('u') 1513 cosmo_lon => cosmo_grid % lonu 1514 cosmo_lat => cosmo_grid % lat 1515 1516 CASE('v') 1517 cosmo_lon => cosmo_grid % lon 1518 cosmo_lat => cosmo_grid % latv 1519 1520 CASE DEFAULT 1521 message = "Averaging grid kind '" // TRIM(avg_grid % kind) // & 1522 "' is not supported. Use 'scalar', 'u', or 'v'." 1523 CALL abort('get_cosmo_averaging_region', message) 1524 1525 END SELECT 1526 1527 1528 ! FIXME: make dlon, dlat parameters of the grid_defintion type 1529 dlon = cosmo_lon(1) - cosmo_lon(0) 1530 dlat = cosmo_lat(1) - cosmo_lat(0) 1531 1532 imin = CEILING( (avg_grid % lon(1) - cosmo_lon(0)) / dlon ) 1533 imax = FLOOR( (avg_grid % lon(2) - cosmo_lon(0)) / dlon ) 1534 1535 jmin = CEILING( (avg_grid % lat(1) - cosmo_lat(0)) / dlat ) 1536 jmax = FLOOR( (avg_grid % lat(2) - cosmo_lat(0)) / dlat ) 1537 1538 message = "Averaging over "// & 1539 TRIM(str(imin)) // " < i < " // TRIM(str(imax)) // & 1540 " and " // & 1541 TRIM(str(jmin)) // " < j < " // TRIM(str(jmax)) 1542 CALL report( 'get_cosmo_averaging_region', message ) 1543 1544 nx = imax - imin + 1 1545 ny = jmax - jmin + 1 1546 avg_grid % n_columns = nx * ny 1547 1548 ALLOCATE( avg_grid % iii(avg_grid % n_columns), & 1549 avg_grid % jjj(avg_grid % n_columns) ) 1550 1551 l = 0 1552 DO j = jmin, jmax 1553 DO i = imin, imax 1554 l = l + 1 1555 avg_grid % iii(l) = i 1556 avg_grid % jjj(l) = j 1557 END DO 1558 END DO 1559 1560 END SUBROUTINE get_cosmo_averaging_region 1293 1561 1294 1562 … … 1718 1986 ) 1719 1987 1720 !potential temperature, surface pressure, including nudging and subsidence 1988 !potential temperature, surface pressure, specific humidity including 1989 !nudging and subsidence, and geostrophic winds ug, vg 1721 1990 io_group_list(3) = init_io_group( & 1722 1991 in_files = flow_files, & 1723 out_vars = [output_var_table(3:8), output_var_table(42:42), & 1724 output_var_table(49:51)], & 1725 in_var_list = (/input_var_table(3), input_var_table(17)/), & 1726 kind = 'temperature' & 1727 ) 1728 1729 !specific humidity including nudging and subsidence 1730 io_group_list(4) = init_io_group( & 1731 in_files = flow_files, & 1732 out_vars = [output_var_table(9:14), output_var_table(52:54)], & 1733 in_var_list = input_var_table(4:4), & 1734 kind = 'scalar' & 1735 ) 1736 1737 !u and v velocity and geostrophic winds ug, vg 1992 out_vars = [output_var_table(56:64), & ! internal averaged density and pressure profiles 1993 output_var_table(3:8), output_var_table(9:14), & 1994 output_var_table(42:42), output_var_table(43:44), & 1995 output_var_table(49:51), output_var_table(52:54)], & 1996 in_var_list = (/input_var_table(3), input_var_table(17), & ! T, P, QV 1997 input_var_table(4) /), & 1998 kind = 'thermodynamics', & 1999 n_output_quantities = 4 & ! P, Theta, Rho, qv 2000 ) 2001 2002 !Moved to therodynamic io_group 2003 !io_group_list(4) = init_io_group( & 2004 ! in_files = flow_files, & 2005 ! out_vars = [output_var_table(9:14), output_var_table(52:54)], & 2006 ! in_var_list = input_var_table(4:4), & 2007 ! kind = 'scalar' & 2008 !) 2009 2010 !u and v velocity, ug anv vg moved to thermodynamic io group. 1738 2011 io_group_list(5) = init_io_group( & 1739 2012 in_files = flow_files, & 1740 out_vars = [output_var_table(15:26), output_var_table(4 3:46)], &2013 out_vars = [output_var_table(15:26), output_var_table(45:46)], & 1741 2014 !out_vars = output_var_table(15:20), & 1742 2015 in_var_list = input_var_table(5:6), & … … 1845 2118 1846 2119 1847 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group) 2120 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind, & 2121 n_output_quantities) RESULT(group) 1848 2122 CHARACTER(LEN=PATH), INTENT(IN) :: in_files(:) 1849 2123 CHARACTER(LEN=*), INTENT(IN) :: kind 1850 2124 TYPE(nc_var), INTENT(IN) :: out_vars(:) 1851 2125 TYPE(nc_var), INTENT(IN) :: in_var_list(:) 2126 INTEGER, OPTIONAL :: n_output_quantities 1852 2127 1853 2128 TYPE(io_group) :: group … … 1855 2130 group % nt = SIZE(in_files) 1856 2131 group % nv = SIZE(out_vars) 2132 group % n_inputs = SIZE(in_var_list) 1857 2133 group % kind = TRIM(kind) 1858 2134 1859 ALLOCATE(group % in_var_list( SIZE(in_var_list) )) 2135 IF ( PRESENT(n_output_quantities) ) THEN 2136 group % n_output_quantities = n_output_quantities 2137 ELSE 2138 group % n_output_quantities = group % n_inputs 2139 END IF 2140 2141 ALLOCATE(group % in_var_list(group % n_inputs)) 1860 2142 ALLOCATE(group % in_files(group % nt)) 1861 2143 ALLOCATE(group % out_vars(group % nv)) … … 1876 2158 SUBROUTINE fini_grids() 1877 2159 1878 CALL report('fini_grids', 'Deallocating grids' )2160 CALL report('fini_grids', 'Deallocating grids', cfg % debug) 1879 2161 1880 2162 DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column) … … 1904 2186 SUBROUTINE setup_variable_tables(ic_mode) 1905 2187 CHARACTER(LEN=*), INTENT(IN) :: ic_mode 2188 INTEGER :: n_invar = 0 !< number of variables in the input variable table 2189 INTEGER :: n_outvar = 0 !< number of variables in the output variable table 1906 2190 TYPE(nc_var), POINTER :: var 1907 2191 … … 1914 2198 1915 2199 n_invar = 17 1916 n_outvar = 552200 n_outvar = 64 1917 2201 ALLOCATE( input_var_table(n_invar) ) 1918 2202 ALLOCATE( output_var_table(n_outvar) ) … … 2052 2336 ) 2053 2337 IF (TRIM(ic_mode) == 'profile') THEN 2054 output_var_table(3) % grid => scalar_profile_grid 2055 output_var_table(3) % intermediate_grid => scalar_profile_intermediate 2338 output_var_table(3) % averaging_grid => averaged_scalar_profile 2056 2339 END IF 2057 2340 … … 2122 2405 units = "kg/kg", & 2123 2406 kind = "init scalar", & 2124 input_id = 1, &2407 input_id = 3, & 2125 2408 output_file = output_file, & 2126 2409 grid = palm_grid, & … … 2129 2412 ) 2130 2413 IF (TRIM(ic_mode) == 'profile') THEN 2131 output_var_table(9) % grid => scalar_profile_grid 2132 output_var_table(9) % intermediate_grid => scalar_profile_intermediate 2414 output_var_table(9) % averaging_grid => averaged_scalar_profile 2133 2415 END IF 2134 2416 … … 2139 2421 units = "kg/kg", & 2140 2422 kind = "left scalar", & 2141 input_id = 1, &2423 input_id = 3, & 2142 2424 output_file = output_file, & 2143 2425 grid = scalars_west_grid, & … … 2151 2433 units = "kg/kg", & 2152 2434 kind = "right scalar", & 2153 input_id = 1, &2435 input_id = 3, & 2154 2436 output_file = output_file, & 2155 2437 grid = scalars_east_grid, & … … 2163 2445 units = "kg/kg", & 2164 2446 kind = "north scalar", & 2165 input_id = 1, &2447 input_id = 3, & 2166 2448 output_file = output_file, & 2167 2449 grid = scalars_north_grid, & … … 2175 2457 units = "kg/kg", & 2176 2458 kind = "south scalar", & 2177 input_id = 1, &2459 input_id = 3, & 2178 2460 output_file = output_file, & 2179 2461 grid = scalars_south_grid, & … … 2187 2469 units = "kg/kg", & 2188 2470 kind = "top scalar", & 2189 input_id = 1, &2471 input_id = 3, & 2190 2472 output_file = output_file, & 2191 2473 grid = scalars_top_grid, & … … 2206 2488 ) 2207 2489 IF (TRIM(ic_mode) == 'profile') THEN 2208 output_var_table(15) % grid => scalar_profile_grid 2209 output_var_table(15) % intermediate_grid => scalar_profile_intermediate 2490 output_var_table(15) % averaging_grid => averaged_scalar_profile 2210 2491 END IF 2211 2492 … … 2283 2564 ) 2284 2565 IF (TRIM(ic_mode) == 'profile') THEN 2285 output_var_table(21) % grid => scalar_profile_grid 2286 output_var_table(21) % intermediate_grid => scalar_profile_intermediate 2566 output_var_table(21) % averaging_grid => averaged_scalar_profile 2287 2567 END IF 2288 2568 … … 2360 2640 ) 2361 2641 IF (TRIM(ic_mode) == 'profile') THEN 2362 output_var_table(27) % grid => w_profile_grid 2363 output_var_table(27) % intermediate_grid => w_profile_intermediate 2642 output_var_table(27) % averaging_grid => averaged_w_profile 2364 2643 END IF 2365 2644 … … 2546 2825 intermediate_grid = palm_intermediate & 2547 2826 ) 2827 output_var_table(42) % averaging_grid => averaged_scalar_profile 2548 2828 2549 2829 output_var_table(43) = init_nc_var( & … … 2552 2832 long_name = "geostrophic wind (u component)", & 2553 2833 units = "m/s", & 2554 kind = " constant scalar profile",&2834 kind = "geostrophic", & 2555 2835 input_id = 1, & 2556 2836 output_file = output_file, & 2557 grid = scalar_profile_grid,&2558 intermediate_grid = scalar_profile_intermediate&2837 grid = averaged_scalar_profile, & 2838 intermediate_grid = averaged_scalar_profile & 2559 2839 ) 2560 2840 … … 2564 2844 long_name = "geostrophic wind (v component)", & 2565 2845 units = "m/s", & 2566 kind = " constant scalar profile",&2846 kind = "geostrophic", & 2567 2847 input_id = 1, & 2568 2848 output_file = output_file, & 2569 grid = scalar_profile_grid,&2570 intermediate_grid = scalar_profile_intermediate&2849 grid = averaged_scalar_profile, & 2850 intermediate_grid = averaged_scalar_profile & 2571 2851 ) 2572 2852 … … 2576 2856 long_name = "wind component in x direction", & 2577 2857 units = "m/s", & 2578 kind = " large-scale scalar forcing",&2858 kind = "geostrophic", & 2579 2859 input_id = 1, & 2580 2860 output_file = output_file, & 2581 grid = scalar_profile_grid,&2582 intermediate_grid = scalar_profile_intermediate&2861 grid = averaged_scalar_profile, & 2862 intermediate_grid = averaged_scalar_profile & 2583 2863 ) 2584 2864 output_var_table(45) % to_be_processed = ls_forcing_variables_required … … 2592 2872 input_id = 1, & 2593 2873 output_file = output_file, & 2594 grid = scalar_profile_grid,&2595 intermediate_grid = scalar_profile_intermediate&2874 grid = averaged_scalar_profile, & 2875 intermediate_grid = averaged_scalar_profile & 2596 2876 ) 2597 2877 output_var_table(46) % to_be_processed = ls_forcing_variables_required … … 2605 2885 input_id = 1, & 2606 2886 output_file = output_file, & 2607 grid = w_profile_grid,&2608 intermediate_grid = w_profile_intermediate&2887 grid = averaged_scalar_profile, & 2888 intermediate_grid = averaged_scalar_profile & 2609 2889 ) 2610 2890 output_var_table(47) % to_be_processed = ls_forcing_variables_required … … 2618 2898 input_id = 1, & 2619 2899 output_file = output_file, & 2620 grid = w_profile_grid,&2621 intermediate_grid = w_profile_intermediate&2900 grid = averaged_w_profile, & 2901 intermediate_grid = averaged_w_profile & 2622 2902 ) 2623 2903 output_var_table(48) % to_be_processed = ls_forcing_variables_required … … 2632 2912 input_id = 1, & 2633 2913 output_file = output_file, & 2634 grid = scalar_profile_grid,&2635 intermediate_grid = scalar_profile_intermediate&2914 grid = averaged_scalar_profile, & 2915 intermediate_grid = averaged_scalar_profile & 2636 2916 ) 2637 2917 output_var_table(49) % to_be_processed = ls_forcing_variables_required … … 2645 2925 input_id = 1, & 2646 2926 output_file = output_file, & 2647 grid = scalar_profile_grid,&2648 intermediate_grid = scalar_profile_intermediate&2927 grid = averaged_scalar_profile, & 2928 intermediate_grid = averaged_scalar_profile & 2649 2929 ) 2650 2930 output_var_table(50) % to_be_processed = ls_forcing_variables_required … … 2658 2938 input_id = 1, & 2659 2939 output_file = output_file, & 2660 grid = scalar_profile_grid,&2661 intermediate_grid = scalar_profile_intermediate&2940 grid = averaged_scalar_profile, & 2941 intermediate_grid = averaged_scalar_profile & 2662 2942 ) 2663 2943 output_var_table(51) % to_be_processed = ls_forcing_variables_required … … 2669 2949 units = "kg/kg/s", & 2670 2950 kind = "large-scale scalar forcing", & 2671 input_id = 1, &2672 output_file = output_file, & 2673 grid = scalar_profile_grid,&2674 intermediate_grid = scalar_profile_intermediate&2951 input_id = 3, & 2952 output_file = output_file, & 2953 grid = averaged_scalar_profile, & 2954 intermediate_grid = averaged_scalar_profile & 2675 2955 ) 2676 2956 output_var_table(52) % to_be_processed = ls_forcing_variables_required … … 2683 2963 units = "kg/kg/s", & 2684 2964 kind = "large-scale scalar forcing", & 2685 input_id = 1, &2686 output_file = output_file, & 2687 grid = scalar_profile_grid,&2688 intermediate_grid = scalar_profile_intermediate&2965 input_id = 3, & 2966 output_file = output_file, & 2967 grid = averaged_scalar_profile, & 2968 intermediate_grid = averaged_scalar_profile & 2689 2969 ) 2690 2970 output_var_table(53) % to_be_processed = ls_forcing_variables_required … … 2696 2976 units = "kg/kg", & 2697 2977 kind = "large-scale scalar forcing", & 2698 input_id = 1, &2699 output_file = output_file, & 2700 grid = scalar_profile_grid,&2701 intermediate_grid = scalar_profile_intermediate&2978 input_id = 3, & 2979 output_file = output_file, & 2980 grid = averaged_scalar_profile, & 2981 intermediate_grid = averaged_scalar_profile & 2702 2982 ) 2703 2983 output_var_table(54) % to_be_processed = ls_forcing_variables_required … … 2711 2991 input_id = 1, & 2712 2992 output_file = output_file, & 2713 grid = scalar_profile_grid,&2714 intermediate_grid = scalar_profile_intermediate&2993 grid = averaged_scalar_profile, & 2994 intermediate_grid = averaged_scalar_profile & 2715 2995 ) 2716 2996 output_var_table(55) % to_be_processed = ls_forcing_variables_required 2717 2997 2998 2999 output_var_table(56) = init_nc_var( & 3000 name = 'internal_density_centre', & 3001 std_name = "", & 3002 long_name = "", & 3003 units = "", & 3004 kind = "internal profile", & 3005 input_id = 4, & 3006 output_file = output_file, & 3007 grid = averaged_scalar_profile, & 3008 intermediate_grid = averaged_scalar_profile & 3009 ) 3010 output_var_table(56) % averaging_grid => averaged_scalar_profile 3011 3012 3013 output_var_table(57) = init_nc_var( & 3014 name = 'internal_density_north', & 3015 std_name = "", & 3016 long_name = "", & 3017 units = "", & 3018 kind = "internal profile", & 3019 input_id = 4, & 3020 output_file = output_file, & 3021 grid = north_averaged_scalar_profile, & 3022 intermediate_grid = north_averaged_scalar_profile & 3023 ) 3024 output_var_table(57) % averaging_grid => north_averaged_scalar_profile 3025 output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set 3026 3027 3028 output_var_table(58) = init_nc_var( & 3029 name = 'internal_density_south', & 3030 std_name = "", & 3031 long_name = "", & 3032 units = "", & 3033 kind = "internal profile", & 3034 input_id = 4, & 3035 output_file = output_file, & 3036 grid = south_averaged_scalar_profile, & 3037 intermediate_grid = south_averaged_scalar_profile & 3038 ) 3039 output_var_table(58) % averaging_grid => south_averaged_scalar_profile 3040 output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set 3041 3042 3043 output_var_table(59) = init_nc_var( & 3044 name = 'internal_density_east', & 3045 std_name = "", & 3046 long_name = "", & 3047 units = "", & 3048 kind = "internal profile", & 3049 input_id = 4, & 3050 output_file = output_file, & 3051 grid = east_averaged_scalar_profile, & 3052 intermediate_grid = east_averaged_scalar_profile & 3053 ) 3054 output_var_table(59) % averaging_grid => east_averaged_scalar_profile 3055 output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set 3056 3057 3058 output_var_table(60) = init_nc_var( & 3059 name = 'internal_density_west', & 3060 std_name = "", & 3061 long_name = "", & 3062 units = "", & 3063 kind = "internal profile", & 3064 input_id = 4, & 3065 output_file = output_file, & 3066 grid = west_averaged_scalar_profile, & 3067 intermediate_grid = west_averaged_scalar_profile & 3068 ) 3069 output_var_table(60) % averaging_grid => west_averaged_scalar_profile 3070 output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set 3071 3072 output_var_table(61) = init_nc_var( & 3073 name = 'internal_pressure_north', & 3074 std_name = "", & 3075 long_name = "", & 3076 units = "", & 3077 kind = "internal profile", & 3078 input_id = 2, & 3079 output_file = output_file, & 3080 grid = north_averaged_scalar_profile, & 3081 intermediate_grid = north_averaged_scalar_profile & 3082 ) 3083 output_var_table(61) % averaging_grid => north_averaged_scalar_profile 3084 output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set 3085 3086 3087 output_var_table(62) = init_nc_var( & 3088 name = 'internal_pressure_south', & 3089 std_name = "", & 3090 long_name = "", & 3091 units = "", & 3092 kind = "internal profile", & 3093 input_id = 2, & 3094 output_file = output_file, & 3095 grid = south_averaged_scalar_profile, & 3096 intermediate_grid = south_averaged_scalar_profile & 3097 ) 3098 output_var_table(62) % averaging_grid => south_averaged_scalar_profile 3099 output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set 3100 3101 3102 output_var_table(63) = init_nc_var( & 3103 name = 'internal_pressure_east', & 3104 std_name = "", & 3105 long_name = "", & 3106 units = "", & 3107 kind = "internal profile", & 3108 input_id = 2, & 3109 output_file = output_file, & 3110 grid = east_averaged_scalar_profile, & 3111 intermediate_grid = east_averaged_scalar_profile & 3112 ) 3113 output_var_table(63) % averaging_grid => east_averaged_scalar_profile 3114 output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set 3115 3116 3117 output_var_table(64) = init_nc_var( & 3118 name = 'internal_pressure_west', & 3119 std_name = "", & 3120 long_name = "", & 3121 units = "", & 3122 kind = "internal profile", & 3123 input_id = 2, & 3124 output_file = output_file, & 3125 grid = west_averaged_scalar_profile, & 3126 intermediate_grid = west_averaged_scalar_profile & 3127 ) 3128 output_var_table(64) % averaging_grid => west_averaged_scalar_profile 3129 output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set 2718 3130 2719 3131 ! Attributes shared among all variables … … 2732 3144 !------------------------------------------------------------------------------! 2733 3145 FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, & 2734 grid, intermediate_grid, output_file, is_profile 2735 )RESULT(var)3146 grid, intermediate_grid, output_file, is_profile) & 3147 RESULT(var) 2736 3148 2737 3149 CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind … … 2772 3184 var % dimvarids(1:3) = output_file % dimvarids_soil 2773 3185 var % to_be_processed = init_variables_required 3186 var % is_internal = .FALSE. 2774 3187 var % task = "interpolate_2d" 2775 3188 … … 2781 3194 var % dimvarids(1:3) = output_file % dimvarids_scl 2782 3195 var % to_be_processed = init_variables_required 3196 var % is_internal = .FALSE. 2783 3197 var % task = "interpolate_3d" 2784 3198 … … 2794 3208 var % dimvarids(3) = output_file % dimvarids_scl(3) 2795 3209 var % to_be_processed = init_variables_required 3210 var % is_internal = .FALSE. 2796 3211 var % task = "interpolate_3d" 2797 3212 … … 2807 3222 var % dimvarids(3) = output_file % dimvarids_scl(3) 2808 3223 var % to_be_processed = init_variables_required 3224 var % is_internal = .FALSE. 2809 3225 var % task = "interpolate_3d" 2810 3226 … … 2820 3236 var % dimvarids(3) = output_file % dimvarids_vel(3) 2821 3237 var % to_be_processed = init_variables_required 3238 var % is_internal = .FALSE. 2822 3239 var % task = "interpolate_3d" 2823 3240 … … 2829 3246 var % dimvarids(1) = output_file % dimvarids_scl(3) !z 2830 3247 var % to_be_processed = init_variables_required 3248 var % is_internal = .FALSE. 2831 3249 var % task = "average profile" 2832 3250 … … 2838 3256 var % dimvarids(1) = output_file % dimvarids_vel(3) !z 2839 3257 var % to_be_processed = init_variables_required 3258 var % is_internal = .FALSE. 2840 3259 var % task = "average profile" 2841 3260 … … 2848 3267 var % dimvarids(1:2) = output_file % dimvarids_soil(1:2) 2849 3268 var % to_be_processed = boundary_variables_required 3269 var % is_internal = .FALSE. 2850 3270 var % task = "interpolate_2d" 2851 3271 … … 2860 3280 var % dimvarids(2) = output_file % dimvarids_scl(3) 2861 3281 var % to_be_processed = boundary_variables_required 3282 var % is_internal = .FALSE. 2862 3283 var % task = "interpolate_3d" 2863 3284 … … 2872 3293 var % dimvarids(2) = output_file % dimvarids_scl(3) 2873 3294 var % to_be_processed = boundary_variables_required 3295 var % is_internal = .FALSE. 2874 3296 var % task = "interpolate_3d" 2875 3297 … … 2884 3306 var % dimvarids(2) = output_file % dimvarids_scl(2) 2885 3307 var % to_be_processed = boundary_variables_required 3308 var % is_internal = .FALSE. 2886 3309 var % task = "interpolate_3d" 2887 3310 … … 2896 3319 var % dimvarids(2) = output_file % dimvarids_scl(3) 2897 3320 var % to_be_processed = boundary_variables_required 3321 var % is_internal = .FALSE. 2898 3322 var % task = "interpolate_3d" 2899 3323 … … 2908 3332 var % dimvarids(2) = output_file % dimvarids_scl(3) 2909 3333 var % to_be_processed = boundary_variables_required 3334 var % is_internal = .FALSE. 2910 3335 var % task = "interpolate_3d" 2911 3336 … … 2920 3345 var % dimvarids(2) = output_file % dimvarids_scl(2) 2921 3346 var % to_be_processed = boundary_variables_required 3347 var % is_internal = .FALSE. 2922 3348 var % task = "interpolate_3d" 2923 3349 … … 2932 3358 var % dimvarids(2) = output_file % dimvarids_scl(3) 2933 3359 var % to_be_processed = boundary_variables_required 3360 var % is_internal = .FALSE. 2934 3361 var % task = "interpolate_3d" 2935 3362 … … 2944 3371 var % dimvarids(2) = output_file % dimvarids_scl(3) 2945 3372 var % to_be_processed = boundary_variables_required 3373 var % is_internal = .FALSE. 2946 3374 var % task = "interpolate_3d" 2947 3375 … … 2956 3384 var % dimvarids(2) = output_file % dimvarids_vel(2) 2957 3385 var % to_be_processed = boundary_variables_required 3386 var % is_internal = .FALSE. 2958 3387 var % task = "interpolate_3d" 2959 3388 … … 2968 3397 var % dimvarids(2) = output_file % dimvarids_vel(3) 2969 3398 var % to_be_processed = boundary_variables_required 3399 var % is_internal = .FALSE. 2970 3400 var % task = "interpolate_3d" 2971 3401 … … 2980 3410 var % dimvarids(2) = output_file % dimvarids_vel(3) 2981 3411 var % to_be_processed = boundary_variables_required 3412 var % is_internal = .FALSE. 2982 3413 var % task = "interpolate_3d" 2983 3414 … … 2988 3419 var % dimvarids(1) = output_file % dimvarid_time 2989 3420 var % to_be_processed = .TRUE. 2990 var % task = "average scalar" 3421 var % is_internal = .FALSE. 3422 var % task = "average profile" 2991 3423 2992 3424 CASE( 'constant scalar profile' ) … … 2998 3430 var % dimvarids(1) = output_file % dimvarids_scl(3) 2999 3431 var % to_be_processed = .TRUE. 3432 var % is_internal = .FALSE. 3000 3433 var % task = "set profile" 3001 3434 … … 3008 3441 var % dimvarids(1) = output_file % dimvarids_scl(3) 3009 3442 var % to_be_processed = ls_forcing_variables_required 3443 var % is_internal = .FALSE. 3010 3444 var % task = "average large-scale profile" 3445 3446 CASE( 'geostrophic' ) 3447 var % lod = -1 3448 var % ndim = 2 3449 var % dimids(2) = output_file % dimid_time !t 3450 var % dimids(1) = output_file % dimids_scl(3) !z 3451 var % dimvarids(2) = output_file % dimvarid_time 3452 var % dimvarids(1) = output_file % dimvarids_scl(3) 3453 var % to_be_processed = .TRUE. 3454 var % is_internal = .FALSE. 3455 var % task = "geostrophic winds" 3011 3456 3012 3457 CASE( 'large-scale w forcing' ) … … 3018 3463 var % dimvarids(1) = output_file % dimvarids_vel(3) 3019 3464 var % to_be_processed = ls_forcing_variables_required 3465 var % is_internal = .FALSE. 3020 3466 var % task = "average large-scale profile" 3467 3468 CASE( 'internal profile' ) 3469 var % lod = -1 3470 var % ndim = 2 3471 var % dimids(2) = output_file % dimid_time !t 3472 var % dimids(1) = output_file % dimids_scl(3) !z 3473 var % dimvarids(2) = output_file % dimvarid_time 3474 var % dimvarids(1) = output_file % dimvarids_scl(3) 3475 var % to_be_processed = .TRUE. 3476 var % is_internal = .TRUE. 3477 var % task = "internal profile" 3021 3478 3022 3479 CASE DEFAULT … … 3031 3488 SUBROUTINE fini_variables() 3032 3489 3033 CALL report('fini_variables', 'Deallocating variable table' )3490 CALL report('fini_variables', 'Deallocating variable table', cfg % debug) 3034 3491 DEALLOCATE( input_var_table ) 3035 3492 … … 3039 3496 SUBROUTINE fini_io_groups() 3040 3497 3041 CALL report('fini_io_groups', 'Deallocating IO groups' )3498 CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug) 3042 3499 DEALLOCATE( io_group_list ) 3043 3500 … … 3047 3504 SUBROUTINE fini_file_lists() 3048 3505 3049 CALL report('fini_file_lists', 'Deallocating file lists' )3506 CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug) 3050 3507 DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files ) 3051 3508 … … 3134 3591 ! rotated velocities 3135 3592 DO k = 1, nz 3136 DO j = 2, ny3137 DO i = 2, nx3593 DO j = 1, ny 3594 DO i = 1, nx 3138 3595 CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k), & 3139 3596 vrot = preprocess_buffer(2) % array(i,j,k), & … … 3171 3628 CALL report('preprocess', message) 3172 3629 3173 CASE( 't emperature' ) ! P and T3630 CASE( 'thermodynamics' ) ! T, P, QV 3174 3631 nx = SIZE(input_buffer(1) % array, 1) 3175 3632 ny = SIZE(input_buffer(1) % array, 2) … … 3192 3649 3193 3650 input_buffer (2) % array(i,j,:) = & 3194 input_buffer (2) % array(i,j,:) + basic_state_pressure(:) 3651 HECTO * input_buffer (2) % array(i,j,:) + & 3652 basic_state_pressure(:) 3195 3653 3196 3654 END DO … … 3204 3662 3205 3663 END IF 3664 ! pressure 3665 input_buffer(2) % is_preprocessed = .TRUE. 3666 3667 ! Copy temperature to last input buffer array 3668 ALLOCATE( & 3669 input_buffer( group % n_output_quantities ) % array (nx, ny, nz) & 3670 ) 3671 input_buffer(group % n_output_quantities) % array(:,:,:) = & 3672 input_buffer(1) % array(:,:,:) 3206 3673 3207 3674 ! Convert absolute to potential temperature 3208 input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) * & 3209 EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) )) 3210 3211 input_buffer(1:2) % is_preprocessed = .TRUE. 3675 CALL potential_temperature( & 3676 t = input_buffer(1) % array(:,:,:), & 3677 p = input_buffer(2) % array(:,:,:), & 3678 p_ref = P_REF, & 3679 r = RD_PALM, & 3680 cp = CP_PALM & 3681 ) 3682 3683 ! potential temperature 3684 input_buffer(1) % is_preprocessed = .TRUE. 3685 3686 ! Convert temperature copy to density 3687 CALL moist_density( & 3688 t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), & 3689 p = input_buffer(2) % array(:,:,:), & 3690 qv = input_buffer(3) % array(:,:,:), & 3691 rd = RD, & 3692 rv = RV & 3693 ) 3694 3695 ! qv 3696 input_buffer(3) % is_preprocessed = .TRUE. 3697 3698 ! density 3699 input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE. 3700 3701 3212 3702 message = "Input buffers for group '" // TRIM(group % kind) // "'"//& 3213 3703 " preprocessed sucessfully." … … 3336 3826 3337 3827 CASE DEFAULT 3338 message = " Bufferkind '" // TRIM(group % kind) // "' is not supported."3828 message = "IO group kind '" // TRIM(group % kind) // "' is not supported." 3339 3829 CALL abort('prerpocess', message) 3340 3830 -
palm/trunk/UTIL/inifor/src/inifor.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Added main loop support for computation of geostrophic winds and surface 29 ! pressure 30 ! Cleaned up terminal output, show some meVssages only with --debug option 31 ! 32 ! 3183 2018-07-27 14:25:55Z suehring 28 33 ! Introduced new PALM grid stretching 29 34 ! Renamend initial-condition mode variable 'mode' to 'ic_mode' … … 56 61 fini_file_lists, preprocess, origin_lon, origin_lat, & 57 62 output_file, io_group_list, output_var_table, & 58 cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,&59 averag e_imin, average_imax, average_jmin, average_jmax60 63 cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3, & 64 averaging_width_ns, averaging_width_ew, phi_n, lambda_n, & 65 lam_centre, phi_centre 61 66 USE io 62 67 USE transform, & 63 ONLY: average_profile, average_2d, interpolate_2d, interpolate_3d 68 ONLY: average_profile, interpolate_2d, interpolate_3d, & 69 geostrophic_winds, extrapolate_density, extrapolate_pressure, & 70 get_surface_pressure 64 71 USE types 65 72 … … 70 77 INTEGER :: iter 71 78 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr 79 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre 80 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_arr 81 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_arr 82 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north 83 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south 84 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east 85 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west 86 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north 87 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south 88 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east 89 REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west 90 REAL(dp), POINTER, DIMENSION(:) :: internal_arr 91 REAL(dp), POINTER, DIMENSION(:) :: ug_vg_arr 72 92 TYPE(nc_var), POINTER :: output_var 73 93 TYPE(io_group), POINTER :: group 74 94 TYPE(container), ALLOCATABLE :: input_buffer(:) 95 LOGICAL, SAVE :: ug_vg_have_been_computed = .FALSE. 96 LOGICAL, SAVE :: debugging_output = .TRUE. 75 97 76 98 !> \mainpage About INIFOR … … 199 221 CALL run_control('time', 'alloc') 200 222 201 202 223 CALL average_profile( & 203 224 input_buffer(output_var % input_id) % array(:,:,:), & 204 output_arr(:,:,:), average_imin, average_imax, & 205 average_jmin, average_jmax, & 206 output_var % intermediate_grid, & 207 output_var % grid) 208 CALL run_control('time', 'comp') 225 output_arr(0,0,:), & 226 output_var % averaging_grid) 227 228 IF ( TRIM(output_var % name) == & 229 'surface_forcing_surface_pressure' ) THEN 230 231 IF ( cfg % p0_is_set ) THEN 232 output_arr(0,0,1) = p0 233 ELSE 234 CALL get_surface_pressure( & 235 output_arr(0,0,:), rho_centre, & 236 output_var % averaging_grid) 237 END IF 238 239 END IF 240 CALL run_control('time', 'comp') 241 242 CASE ( 'internal profile' ) 243 244 message = "Averaging of internal profile for variable '" //& 245 TRIM(output_var % name) // "' is not supported." 246 247 SELECT CASE (TRIM(output_var % name)) 248 249 CASE('internal_density_centre') 250 ALLOCATE( rho_centre( 1:output_var % grid % nz ) ) 251 internal_arr => rho_centre 252 253 CASE('internal_density_north') 254 ALLOCATE( rho_north( 1:output_var % grid % nz ) ) 255 internal_arr => rho_north 256 257 CASE('internal_density_south') 258 ALLOCATE( rho_south( 1:output_var % grid % nz ) ) 259 internal_arr => rho_south 260 261 CASE('internal_density_east') 262 ALLOCATE( rho_east( 1:output_var % grid % nz) ) 263 internal_arr => rho_east 264 265 CASE('internal_density_west') 266 ALLOCATE( rho_west( 1:output_var % grid % nz ) ) 267 internal_arr => rho_west 268 269 CASE('internal_pressure_north') 270 ALLOCATE( p_north( 1:output_var % grid % nz ) ) 271 internal_arr => p_north 272 273 CASE('internal_pressure_south') 274 ALLOCATE( p_south( 1:output_var % grid % nz ) ) 275 internal_arr => p_south 276 277 CASE('internal_pressure_east') 278 ALLOCATE( p_east( 1:output_var % grid % nz) ) 279 internal_arr => p_east 280 281 CASE('internal_pressure_west') 282 ALLOCATE( p_west( 1:output_var % grid % nz ) ) 283 internal_arr => p_west 284 285 CASE DEFAULT 286 CALL abort('main loop', message) 287 288 END SELECT 289 CALL run_control('time', 'alloc') 290 291 292 CALL average_profile( & 293 input_buffer(output_var % input_id) % array(:,:,:),& 294 internal_arr(:), & 295 output_var % averaging_grid) 296 297 SELECT CASE (TRIM(output_var % name)) 298 299 CASE('internal_density_centre', & 300 'internal_density_north', & 301 'internal_density_south', & 302 'internal_density_east', & 303 'internal_density_west') 304 CALL extrapolate_density(internal_arr, & 305 output_var % averaging_grid) 306 307 CASE('internal_pressure_north') 308 CALL extrapolate_pressure(internal_arr, rho_north, & 309 output_var % averaging_grid) 310 311 CASE('internal_pressure_south') 312 CALL extrapolate_pressure(internal_arr, rho_south, & 313 output_var % averaging_grid) 314 315 CASE('internal_pressure_east') 316 CALL extrapolate_pressure(internal_arr, rho_east, & 317 output_var % averaging_grid) 318 319 CASE('internal_pressure_west') 320 CALL extrapolate_pressure(internal_arr, rho_west, & 321 output_var % averaging_grid) 322 323 CASE DEFAULT 324 CALL abort('main loop', message) 325 326 END SELECT 327 328 IF (.TRUE.) THEN 329 ALLOCATE( output_arr(1,1,1:output_var % grid % nz) ) 330 output_arr(1,1,:) = internal_arr(:) 331 END IF 332 CALL run_control('time', 'comp') 333 334 335 ! This case gets called twice, the first time for ug, the 336 ! second time for vg. We compute ug and vg at the first call 337 ! and keep vg (and ug for that matter) around for the second 338 ! call. 339 CASE ( 'geostrophic winds' ) 340 341 IF (.NOT. ug_vg_have_been_computed ) THEN 342 ALLOCATE( ug_arr(output_var % grid % nz) ) 343 ALLOCATE( vg_arr(output_var % grid % nz) ) 344 345 IF ( cfg % ug_is_set ) THEN 346 ug_arr = cfg % ug 347 vg_arr = cfg % vg 348 ELSE 349 CALL geostrophic_winds( p_north, p_south, p_east, & 350 p_west, rho_centre, f3, & 351 averaging_width_ew, & 352 averaging_width_ns, & 353 phi_n, lambda_n, phi_centre, lam_centre, & 354 ug_arr, vg_arr ) 355 END IF 356 357 ug_vg_have_been_computed = .TRUE. 358 359 END IF 360 361 ! Prepare output of geostrophic winds 362 SELECT CASE(TRIM(output_var % name)) 363 CASE ('ls_forcing_ug') 364 ug_vg_arr => ug_arr 365 CASE ('ls_forcing_vg') 366 ug_vg_arr => vg_arr 367 END SELECT 368 369 ALLOCATE( output_arr(1,1,output_var % grid % nz) ) 370 output_arr(1,1,:) = ug_vg_arr(:) 209 371 210 372 CASE ( 'average scalar' ) … … 222 384 SELECT CASE (TRIM(output_var % name)) 223 385 224 CASE('ls_forcing_ug')225 output_arr(1, 1, :) = ug226 227 CASE('ls_forcing_vg')228 output_arr(1, 1, :) = vg386 !CASE('ls_forcing_ug') 387 ! output_arr(1, 1, :) = ug 388 389 !CASE('ls_forcing_vg') 390 ! output_arr(1, 1, :) = vg 229 391 230 392 CASE('nudging_tau') … … 256 418 !- Section 2.3: Write current time step of current variable 257 419 !------------------------------------------------------------------------------ 258 message = "Writing variable '" // TRIM(output_var%name) // "'." 259 CALL report('main loop', message) 260 CALL update_output(output_var, output_arr, iter, output_file) 420 IF (.NOT. output_var % is_internal .OR. debugging_output) THEN 421 message = "Writing variable '" // TRIM(output_var%name) // "'." 422 CALL report('main loop', message) 423 CALL update_output(output_var, output_arr, iter, output_file) 261 424 CALL run_control('time', 'write') 262 263 DEALLOCATE(output_arr) 425 END IF 426 427 IF (ALLOCATED(output_arr)) DEALLOCATE(output_arr) 264 428 CALL run_control('time', 'alloc') 265 429 … … 267 431 268 432 END DO ! ouput variables 433 434 IF ( group % kind == 'thermodynamics' ) THEN 435 DEALLOCATE( rho_centre ) 436 IF ( .NOT. cfg % ug_is_set ) THEN 437 DEALLOCATE( rho_north ) 438 DEALLOCATE( rho_south ) 439 DEALLOCATE( rho_east ) 440 DEALLOCATE( rho_west ) 441 DEALLOCATE( p_north ) 442 DEALLOCATE( p_south ) 443 DEALLOCATE( p_east ) 444 DEALLOCATE( p_west ) 445 END IF 446 END IF 269 447 270 448 IF ( group % kind == 'running average' .OR. & … … 273 451 ! accumulated COSMO-DE quantities (precipitation). 274 452 ELSE 275 CALL report('main loop', 'Deallocating input buffer' )453 CALL report('main loop', 'Deallocating input buffer', cfg % debug) 276 454 DEALLOCATE(input_buffer) 277 455 END IF … … 281 459 282 460 IF (ALLOCATED(input_buffer)) THEN 283 CALL report('main loop', 'Deallocating input buffer' )461 CALL report('main loop', 'Deallocating input buffer', cfg % debug) 284 462 DEALLOCATE(input_buffer) 285 463 END IF … … 294 472 END IF 295 473 296 CALL report('main loop', message )474 CALL report('main loop', message, cfg % debug) 297 475 298 476 END IF ! IO group % to_be_processed -
palm/trunk/UTIL/inifor/src/io.f90
r3262 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Removed unnecessary check for output file 28 ! Added command-line options for configuring the computation of geostrophic 29 ! winds (--averagin-mode, --averaging-angle) 30 ! Added command-line option --input-prefix for setting input file prefixes all 31 ! at once 32 ! Added --debug option for more verbose terminal output 33 ! Added option-specific *_is_set LOGICALs to indicate invocation from the 34 ! command-line 35 ! Improved error messages in case of empty file-name strings 36 ! Improved routine naming 29 37 ! 30 38 ! 3183 2018-07-27 14:25:55Z suehring … … 148 156 INTEGER :: arg_count, i 149 157 158 cfg % p0_is_set = .FALSE. 159 cfg % ug_is_set = .FALSE. 160 cfg % vg_is_set = .FALSE. 161 cfg % flow_prefix_is_set = .FALSE. 162 cfg % input_prefix_is_set = .FALSE. 163 cfg % radiation_prefix_is_set = .FALSE. 164 cfg % soil_prefix_is_set = .FALSE. 165 cfg % soilmoisture_prefix_is_set = .FALSE. 166 150 167 arg_count = COMMAND_ARGUMENT_COUNT() 151 168 IF (arg_count .GT. 0) THEN … … 170 187 SELECT CASE( TRIM(option) ) 171 188 189 CASE( '--averaging-mode' ) 190 CALL get_option_argument( i, arg ) 191 cfg % averaging_mode = TRIM(arg) 192 172 193 CASE( '-date', '-d', '--date' ) 173 194 CALL get_option_argument( i, arg ) 174 195 cfg % start_date = TRIM(arg) 196 197 CASE( '--debug' ) 198 cfg % debug = .TRUE. 175 199 176 200 ! Elevation of the PALM-4U domain above sea level … … 181 205 ! surface pressure, at z0 182 206 CASE( '-p0', '-r', '--surface-pressure' ) 207 cfg % p0_is_set = .TRUE. 183 208 CALL get_option_argument( i, arg ) 184 209 READ(arg, *) cfg % p0 … … 186 211 ! geostrophic wind in x direction 187 212 CASE( '-ug', '-u', '--geostrophic-u' ) 213 cfg % ug_is_set = .TRUE. 188 214 CALL get_option_argument( i, arg ) 189 215 READ(arg, *) cfg % ug … … 191 217 ! geostrophic wind in y direction 192 218 CASE( '-vg', '-v', '--geostrophic-v' ) 219 cfg % vg_is_set = .TRUE. 193 220 CALL get_option_argument( i, arg ) 194 221 READ(arg, *) cfg % vg … … 208 235 CASE( '-hhl', '-l', '--hhl-file' ) 209 236 CALL get_option_argument( i, arg ) 210 cfg % hhl_file = TRIM(arg) 237 cfg % hhl_file = TRIM(arg) 238 239 CASE( '--input-prefix') 240 CALL get_option_argument( i, arg ) 241 cfg % input_prefix = TRIM(arg) 242 cfg % input_prefix_is_set = .TRUE. 243 244 CASE( '-a', '--averaging-angle' ) 245 CALL get_option_argument( i, arg ) 246 READ(arg, *) cfg % averaging_angle 211 247 212 248 CASE( '-static', '-t', '--static-driver' ) 213 249 CALL get_option_argument( i, arg ) 214 250 cfg % static_driver_file = TRIM(arg) 215 251 216 252 CASE( '-soil', '-s', '--soil-file') 217 253 CALL get_option_argument( i, arg ) 218 254 cfg % soiltyp_file = TRIM(arg) 219 255 220 256 CASE( '--flow-prefix') 221 257 CALL get_option_argument( i, arg ) 222 cfg % flow_prefix = TRIM(arg) 223 258 cfg % flow_prefix = TRIM(arg) 259 cfg % flow_prefix_is_set = .TRUE. 260 224 261 CASE( '--radiation-prefix') 225 262 CALL get_option_argument( i, arg ) 226 cfg % radiation_prefix = TRIM(arg) 227 263 cfg % radiation_prefix = TRIM(arg) 264 cfg % radiation_prefix_is_set = .TRUE. 265 228 266 CASE( '--soil-prefix') 229 267 CALL get_option_argument( i, arg ) 230 cfg % soil_prefix = TRIM(arg) 231 268 cfg % soil_prefix = TRIM(arg) 269 cfg % soil_prefix_is_set = .TRUE. 270 232 271 CASE( '--soilmoisture-prefix') 233 272 CALL get_option_argument( i, arg ) 234 cfg % soilmoisture_prefix = TRIM(arg) 273 cfg % soilmoisture_prefix = TRIM(arg) 274 cfg % soilmoisture_prefix_is_set = .TRUE. 235 275 236 276 CASE( '-o', '--output' ) … … 297 337 298 338 all_files_present = .TRUE. 299 all_files_present = all_files_present .AND. file_present(cfg % hhl_file )300 all_files_present = all_files_present .AND. file_present(cfg % namelist_file )301 all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file )339 all_files_present = all_files_present .AND. file_present(cfg % hhl_file, 'HHL') 340 all_files_present = all_files_present .AND. file_present(cfg % namelist_file, 'NAMELIST') 341 all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP') 302 342 303 343 ! Only check optional static driver file name, if it has been given. 304 344 IF (TRIM(cfg % static_driver_file) .NE. '') THEN 305 all_files_present = all_files_present .AND. file_present(cfg % static_driver_file )345 all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver') 306 346 END IF 307 347 … … 335 375 END SELECT 336 376 377 SELECT CASE( TRIM(cfg % averaging_mode) ) 378 CASE( 'level', 'height') 379 CASE DEFAULT 380 message = "Averaging mode '" // TRIM(cfg % averaging_mode) //& 381 "' is not supported. " //& 382 "Please select either 'height' or 'level', " //& 383 "or omit the --averaging-mode option entirely, which corresponds "//& 384 "to the latter." 385 CALL abort( 'validate_config', message ) 386 END SELECT 387 388 IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set ) THEN 389 message = "You specified only one component of the geostrophic " // & 390 "wind. Please specify either both or none." 391 CALL abort( 'validate_config', message ) 392 END IF 337 393 338 394 END SUBROUTINE validate_config 339 395 340 396 341 LOGICAL FUNCTION file_present(filename )397 LOGICAL FUNCTION file_present(filename, kind) 342 398 CHARACTER(LEN=PATH), INTENT(IN) :: filename 343 344 INQUIRE(FILE=filename, EXIST=file_present) 345 346 IF (.NOT. file_present) THEN 347 message = "The given file '" // "' does not exist." 399 CHARACTER(LEN=*), INTENT(IN) :: kind 400 401 IF (LEN(TRIM(filename))==0) THEN 402 403 file_present = .FALSE. 404 message = "No name was given for the " // TRIM(kind) // " file." 348 405 CALL report('file_present', message) 406 407 ELSE 408 409 INQUIRE(FILE=filename, EXIST=file_present) 410 411 IF (.NOT. file_present) THEN 412 message = "The given " // TRIM(kind) // " file '" // & 413 TRIM(filename) // "' does not exist." 414 CALL report('file_present', message) 415 END IF 416 349 417 END IF 350 418 … … 388 456 #endif 389 457 390 ! 391 !------------------------------------------------------------------------------ 392 !- Section 1: Write global NetCDF attributes 458 !------------------------------------------------------------------------------ 459 !- Section 1: Define NetCDF dimensions and coordinates 460 !------------------------------------------------------------------------------ 461 nt = SIZE(output_file % time) 462 nx = palm_grid % nx 463 ny = palm_grid % ny 464 nz = palm_grid % nz 465 z0 = palm_grid % z0 466 467 468 ! 469 !------------------------------------------------------------------------------ 470 !- Section 2: Write global NetCDF attributes 393 471 !------------------------------------------------------------------------------ 394 472 CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string) … … 406 484 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)')))) 407 485 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)')))) 486 ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does 487 ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data 488 ! FIXME: Standard v1.9., origin_z) 489 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)')))) 408 490 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION))) 409 491 CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--')) 410 492 411 493 ! 412 !------------------------------------------------------------------------------413 !- Section 2: Define NetCDF dimensions and coordinates414 !------------------------------------------------------------------------------415 nt = SIZE(output_file % time)416 nx = palm_grid % nx417 ny = palm_grid % ny418 nz = palm_grid % nz419 z0 = palm_grid % z0420 421 494 ! 422 495 !------------------------------------------------------------------------------ … … 618 691 ELSE 619 692 620 nbuffers = SIZE( group % in_var_list ) 693 ! Allocate one input buffer per input_variable. If more quantities 694 ! have to be computed than input variables exist in this group, 695 ! allocate more buffers. For instance, in the thermodynamics group, 696 ! there are three input variabels (PP, T, Qv) and four quantities 697 ! necessart (P, Theta, Rho, qv) for the corresponding output fields 698 ! (p0, Theta, qv, ug, and vg) 699 nbuffers = MAX( group % n_inputs, group % n_output_quantities ) 621 700 ALLOCATE( buffer(nbuffers) ) 622 701 CALL run_control('time', 'alloc') 623 702 624 DO ivar = 1, nbuffers 703 ! Read in all input variables, leave extra buffers-if any-untouched. 704 DO ivar = 1, group % n_inputs 625 705 626 706 input_var => group % in_var_list(ivar) … … 628 708 ! Check wheather P or PP is present in input file 629 709 IF (input_var % name == 'P') THEN 630 input_var % name = TRIM( get_pressure_var (input_file) )710 input_var % name = TRIM( get_pressure_varname(input_file) ) 631 711 CALL run_control('time', 'read') 632 712 END IF … … 668 748 !> perturbation, 'PP', and returns the appropriate string. 669 749 !------------------------------------------------------------------------------! 670 CHARACTER(LEN=2) FUNCTION get_pressure_var (input_file) RESULT(var)750 CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var) 671 751 CHARACTER(LEN=*) :: input_file 672 752 INTEGER :: ncid, varid … … 692 772 CALL check(nf90_close(ncid)) 693 773 694 END FUNCTION get_pressure_var 774 END FUNCTION get_pressure_varname 695 775 696 776 … … 805 885 start=start(1:ndim+1) ) ) 806 886 807 CASE ( 'constant scalar profile' )887 CASE ( 'constant scalar profile', 'geostrophic', 'internal profile' ) 808 888 809 889 CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), & -
palm/trunk/UTIL/inifor/src/transform.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Switched addressing of averaging regions from index bounds to list of columns 29 ! Added routines for the computation of geostrophic winds including: 30 ! - the downward extrapolation of density (linear) and 31 ! - pressure (hydrostatic equation) and 32 ! - rotation of geostrophic wind components to PALM frame of reference 33 ! 34 ! 3183 2018-07-27 14:25:55Z suehring 28 35 ! Introduced new PALM grid stretching 29 36 ! Removed unnecessary subroutine parameters … … 51 58 USE control 52 59 USE defs, & 53 ONLY: TO_DEGREES, TO_RADIANS, PI, dp60 ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp 54 61 USE types 55 62 USE util, & … … 176 183 177 184 178 SUBROUTINE average_2d(in_arr, out_arr, i min, imax, jmin, jmax)179 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:)180 REAL(dp), INTENT(OUT) :: out_arr(0:)181 INTEGER, INTENT(IN) :: imin, imax, jmin, jmax182 183 INTEGER :: i, j, k 185 SUBROUTINE average_2d(in_arr, out_arr, ii, jj) 186 REAL(dp), INTENT(IN) :: in_arr(0:,0:,0:) 187 REAL(dp), INTENT(OUT) :: out_arr(0:) 188 INTEGER, INTENT(IN), DIMENSION(:) :: ii, jj 189 190 INTEGER :: i, j, k, l 184 191 REAL(dp) :: ni 185 186 IF (imin < 0) CALL abort('average_2d', "imin < 0.") 187 IF (jmin < 0) CALL abort('average_2d', "jmin < 0.") 188 IF (imax > UBOUND(in_arr, 1)) CALL abort('average_2d', "imax out of i bound.") 189 IF (jmax > UBOUND(in_arr, 2)) CALL abort('average_2d', "jmax out of j bound.") 192 193 IF (SIZE(ii) .NE. SIZE(jj)) THEN 194 message = "Length of 'ii' and 'jj' index lists do not match." // & 195 NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " // & 196 NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "." 197 CALL abort('average_2d', message) 198 END IF 190 199 191 200 DO k = 0, UBOUND(out_arr, 1) 192 201 193 202 out_arr(k) = 0.0_dp 194 DO j = jmin, jmax 195 DO i = imin, imax 196 out_arr(k) = out_arr(k) + in_arr(i, j, k) 203 DO l = 1, UBOUND(ii, 1) 204 i = ii(l) 205 j = jj(l) 206 out_arr(k) = out_arr(k) +& 207 in_arr(i, j, k) 197 208 END DO 198 END DO 199 200 END DO 201 202 ! devide by number of grid points 203 ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) ) 209 210 END DO 211 212 ni = 1.0_dp / SIZE(ii) 204 213 out_arr(:) = out_arr(:) * ni 205 214 … … 212 221 REAL(dp), DIMENSION(:,:,:), INTENT(OUT) :: palm_array 213 222 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array 214 INTEGER :: nx, ny, n z223 INTEGER :: nx, ny, nlev 215 224 216 225 nx = palm_intermediate % nx 217 226 ny = palm_intermediate % ny 218 n z= palm_intermediate % nz ! nlev227 nlev = palm_intermediate % nz ! nlev 219 228 220 229 ! Interpolate from COSMO-DE to intermediate grid. Allocating with one 221 230 ! less point in the vertical, since scalars like T have 50 instead of 51 222 231 ! points in COSMO-DE. 223 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:n z-1)) !232 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) ! 224 233 225 234 CALL interpolate_2d(source_array, intermediate_array, palm_intermediate) … … 234 243 235 244 236 SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,& 237 palm_intermediate, palm_grid) 238 TYPE(grid_definition), INTENT(IN) :: palm_intermediate, palm_grid 245 SUBROUTINE average_profile(source_array, profile_array, avg_grid) 246 TYPE(grid_definition), INTENT(IN) :: avg_grid 239 247 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: source_array 240 INTEGER, INTENT(IN) :: imin, imax, jmin, jmax 241 REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) :: profile_array 242 REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: intermediate_array 243 INTEGER :: nx, ny, nz 244 245 nx = palm_intermediate % nx 246 ny = palm_intermediate % ny 247 nz = palm_intermediate % nz 248 ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) 249 intermediate_array(:,:,:) = 0.0_dp 250 251 ! average input array to intermediate profile 252 CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax) 253 254 ! vertically interpolate to ouput array 255 CALL interpolate_1d(intermediate_array, profile_array, palm_grid) 256 257 DEALLOCATE(intermediate_array) 248 REAL(dp), DIMENSION(:), INTENT(OUT) :: profile_array 249 250 INTEGER :: i_source, j_source, k_profile, k_source, l, m 251 252 REAL :: ni_columns 253 254 profile_array(:) = 0.0_dp 255 256 DO l = 1, avg_grid % n_columns 257 i_source = avg_grid % iii(l) 258 j_source = avg_grid % jjj(l) 259 260 DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels 261 262 DO m = 1, 2 ! vertical interpolation neighbours 263 264 k_source = avg_grid % kkk(l, k_profile, m) 265 266 profile_array(k_profile) = profile_array(k_profile) & 267 + avg_grid % w(l, k_profile, m) & 268 * source_array(i_source, j_source, k_source) 269 270 END DO ! m, vertical interpolation neighbours 271 272 END DO ! k_profile, PALM levels 273 274 END DO ! l, horizontal neighbours 275 276 ni_columns = 1.0_dp / avg_grid % n_columns 277 profile_array(:) = profile_array(:) * ni_columns 278 279 ! Extrapolate constant to the bottom 280 profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min) 258 281 259 282 END SUBROUTINE average_profile 260 283 284 285 SUBROUTINE extrapolate_density(rho, avg_grid) 286 REAL(dp), DIMENSION(:), INTENT(INOUT) :: rho 287 TYPE(grid_definition), INTENT(IN) :: avg_grid 288 289 REAL(dp) :: drhodz, dz, zk, rhok 290 INTEGER :: k_min 291 292 k_min = avg_grid % k_min 293 zk = avg_grid % z(k_min) 294 rhok = rho(k_min) 295 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 296 drhodz = (rho(k_min + 1) - rho(k_min)) / dz 297 298 rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk) 299 300 END SUBROUTINE extrapolate_density 301 302 303 SUBROUTINE extrapolate_pressure(p, rho, avg_grid) 304 REAL(dp), DIMENSION(:), INTENT(IN) :: rho 305 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p 306 TYPE(grid_definition), INTENT(IN) :: avg_grid 307 308 REAL(dp) :: drhodz, dz, zk, rhok 309 INTEGER :: k, k_min 310 311 k_min = avg_grid % k_min 312 zk = avg_grid % z(k_min) 313 rhok = rho(k_min) 314 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 315 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz 316 317 DO k = 1, k_min-1 318 p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 319 avg_grid % z(k), G) 320 END DO 321 322 END SUBROUTINE extrapolate_pressure 323 324 325 !------------------------------------------------------------------------------! 326 ! Description: 327 ! ------------ 328 !> Takes the averaged pressure profile <p> and sets the lowest entry to the 329 !> extrapolated pressure at the surface. 330 !------------------------------------------------------------------------------! 331 SUBROUTINE get_surface_pressure(p, rho, avg_grid) 332 REAL(dp), DIMENSION(:), INTENT(IN) :: rho 333 REAL(dp), DIMENSION(:), INTENT(INOUT) :: p 334 TYPE(grid_definition), INTENT(IN) :: avg_grid 335 336 REAL(dp) :: drhodz, dz, zk, rhok 337 INTEGER :: k, k_min 338 339 k_min = avg_grid % k_min 340 zk = avg_grid % z(k_min) 341 rhok = rho(k_min) 342 dz = avg_grid % z(k_min + 1) - avg_grid % z(k_min) 343 drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz 344 345 p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz, & 346 0.0, G) 347 348 END SUBROUTINE get_surface_pressure 349 350 351 FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g) RESULT(p) 352 353 REAL(dp), INTENT(IN) :: pk, zk, rhok, drhodz, g, z 354 REAL(dp) :: p 355 356 p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) ) 357 358 END FUNCTION constant_density_pressure 359 360 !-----------------------------------------------------------------------------! 361 ! Description: 362 ! ----------- 363 !> This routine computes profiles of the two geostrophic wind components ug and 364 !> vg. 365 !-----------------------------------------------------------------------------! 366 SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3, & 367 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg) 368 369 REAL(dp), DIMENSION(:), INTENT(IN) :: p_north, p_south, p_east, p_west, & 370 rho 371 REAL(dp), INTENT(IN) :: f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g 372 REAL(dp), DIMENSION(:), INTENT(OUT) :: ug, vg 373 REAL(dp) :: facx, facy 374 375 facx = 1.0_dp / (Lx * f3) 376 facy = 1.0_dp / (Ly * f3) 377 ug(:) = - facy / rho(:) * (p_north(:) - p_south(:)) 378 vg(:) = facx / rho(:) * (p_east(:) - p_west(:)) 379 380 CALL rotate_vector_field( & 381 ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)& 382 ) 383 384 END SUBROUTINE geostrophic_winds 261 385 262 386 … … 419 543 420 544 545 SUBROUTINE rotate_vector_field(x, y, angle) 546 REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y !< x and y coodrinate in arbitrary units 547 REAL(dp), INTENT(IN) :: angle !< rotation angle [deg] 548 549 INTEGER :: i 550 REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2) 551 552 sine = SIN(angle * TO_RADIANS) 553 cosine = COS(angle * TO_RADIANS) 554 ! RESAHPE() fills columns first, so the rotation matrix becomes 555 ! 556 ! rotation = [ cosine -sine ] 557 ! [ sine cosine ] 558 rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) ) 559 560 DO i = LBOUND(x, 1), UBOUND(x, 1) 561 562 v_rot(:) = MATMUL(rotation, (/x(i), y(i)/)) 563 564 x(i) = v_rot(1) 565 y(i) = v_rot(2) 566 567 END DO 568 569 END SUBROUTINE rotate_vector_field 570 571 572 573 !------------------------------------------------------------------------------! 574 ! Description: 575 ! ------------ 576 !> This routine computes the local meridian convergence between a rotated-pole 577 !> and a geographical system using the Eq. (6) given in the DWD manual 578 !> 579 !> Baldauf et al. (2018), Beschreibung des operationelle KuÌrzestfrist- 580 !> vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die 581 !> Datenbanken des DWD. 582 !> https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2 583 !> 584 FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g) & 585 RESULT(delta) 586 587 REAL(dp), INTENT(IN) :: phi_n, lam_n, phi_g, lam_g 588 REAL(dp) :: delta 589 590 delta = atan2( COS(phi_n) * SIN(lam_n - lam_g), & 591 COS(phi_g) * SIN(phi_n) - & 592 SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) ) 593 594 END FUNCTION meridian_convergence_rotated 421 595 422 596 !------------------------------------------------------------------------------! … … 512 686 513 687 514 SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate) 688 SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid, & 689 palm_intermediate ) 515 690 TYPE(grid_definition), INTENT(INOUT) :: palm_grid 516 691 TYPE(grid_definition), INTENT(IN) :: palm_intermediate … … 528 703 529 704 ! in each column of the fine grid, find vertical neighbours of every cell 705 DO j = 0, ny 530 706 DO i = 0, nx 531 DO j = 0, ny532 707 533 708 k_intermediate = 0 … … 536 711 column_top = palm_intermediate % h(i,j,nlev) 537 712 538 ! scan through palm_grid column untiland set neighbour indices in713 ! scan through palm_grid column and set neighbour indices in 539 714 ! case current_height is either below column_base, in the current 540 715 ! cell, or above column_top. Keep increasing current cell index until … … 580 755 h_top = palm_intermediate % h(i,j,k_intermediate+1) 581 756 h_bottom = palm_intermediate % h(i,j,k_intermediate) 757 point_is_in_current_cell = ( & 758 current_height >= h_bottom .AND. & 759 current_height < h_top & 760 ) 761 END DO 762 763 IF (k_intermediate > nlev-1) THEN 764 message = "Index " // TRIM(str(k_intermediate)) // & 765 " is above intermediate grid range." 766 CALL abort('find_vertical_neighbours', message) 767 END IF 768 769 palm_grid % kk(i,j,k,1) = k_intermediate 770 palm_grid % kk(i,j,k,2) = k_intermediate + 1 771 772 ! copmute vertical weights 773 weight = (h_top - current_height) / (h_top - h_bottom) 774 palm_grid % w_verti(i,j,k,1) = weight 775 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight 776 END IF 777 778 END DO 779 780 END DO 781 END DO 782 783 END SUBROUTINE find_vertical_neighbours_and_weights_interp 784 785 786 SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid ) 787 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 788 789 INTEGER :: i, j, k_palm, k_intermediate, l, nlev 790 LOGICAL :: point_is_below_grid, point_is_above_grid, & 791 point_is_in_current_cell 792 REAL(dp) :: current_height, column_base, column_top, h_top, h_bottom, & 793 weight 794 795 796 avg_grid % k_min = LBOUND(avg_grid % z, 1) 797 798 nlev = SIZE(avg_grid % cosmo_h, 3) 799 800 ! in each column of the fine grid, find vertical neighbours of every cell 801 DO l = 1, avg_grid % n_columns 802 803 i = avg_grid % iii(l) 804 j = avg_grid % jjj(l) 805 806 column_base = avg_grid % cosmo_h(i,j,1) 807 column_top = avg_grid % cosmo_h(i,j,nlev) 808 809 ! scan through avg_grid column until and set neighbour indices in 810 ! case current_height is either below column_base, in the current 811 ! cell, or above column_top. Keep increasing current cell index until 812 ! the current cell overlaps with the current_height. 813 k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based. 814 DO k_palm = 1, avg_grid % nz 815 816 ! Memorize the top and bottom boundaries of the coarse cell and the 817 ! current height within it 818 current_height = avg_grid % z(k_palm) + avg_grid % z0 819 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) 820 h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) 821 822 point_is_above_grid = (current_height > column_top) !22000m, very unlikely 823 point_is_below_grid = (current_height < column_base) 824 825 point_is_in_current_cell = ( & 826 current_height >= h_bottom .AND. & 827 current_height < h_top & 828 ) 829 830 ! set default weights 831 avg_grid % w(l,k_palm,1:2) = 0.0_dp 832 833 IF (point_is_above_grid) THEN 834 835 avg_grid % kkk(l,k_palm,1:2) = nlev 836 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 837 838 message = "PALM-4U grid extends above COSMO-DE model top." 839 CALL abort('find_vertical_neighbours_and_weights_average', message) 840 841 ELSE IF (point_is_below_grid) THEN 842 843 avg_grid % kkk(l,k_palm,1:2) = 0 844 avg_grid % w(l,k_palm,1:2) = - 2.0_dp 845 avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min) 846 ELSE 847 ! cycle through intermediate levels until current 848 ! intermediate-grid cell overlaps with current_height 849 DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1) 850 k_intermediate = k_intermediate + 1 851 852 h_top = avg_grid % cosmo_h(i,j,k_intermediate+1) 853 h_bottom = avg_grid % cosmo_h(i,j,k_intermediate) 582 854 point_is_in_current_cell = ( & 583 855 current_height >= h_bottom .AND. & … … 594 866 END IF 595 867 596 palm_grid % kk(i,j,k,1) = k_intermediate597 palm_grid % kk(i,j,k,2) = k_intermediate + 1868 avg_grid % kkk(l,k_palm,1) = k_intermediate 869 avg_grid % kkk(l,k_palm,2) = k_intermediate + 1 598 870 599 871 ! copmute vertical weights 600 872 weight = (h_top - current_height) / (h_top - h_bottom) 601 palm_grid % w_verti(i,j,k,1) = weight602 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight873 avg_grid % w(l,k_palm,1) = weight 874 avg_grid % w(l,k_palm,2) = 1.0_dp - weight 603 875 END IF 604 876 605 END DO 606 607 END DO 608 END DO 609 610 END SUBROUTINE find_vertical_neighbours_and_weights 877 END DO ! k, PALM levels 878 END DO ! l, averaging columns 879 880 END SUBROUTINE find_vertical_neighbours_and_weights_average 611 881 612 882 !------------------------------------------------------------------------------! -
palm/trunk/UTIL/inifor/src/types.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Added *_is_set LOGICALs to inifor_config type to indicate option invocation 29 ! from the command-line 30 ! Added 1D index vertical weights lists to support addressing averaging regions 31 ! by list of columns instead of index bounds 32 ! 33 ! 34 ! 3183 2018-07-27 14:25:55Z suehring 28 35 ! Introduced new PALM grid stretching: 29 36 ! - Converted vertical grid_definition coordinte variables to pointers … … 66 73 67 74 CHARACTER(LEN=SNAME) :: flow_prefix !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses 75 CHARACTER(LEN=SNAME) :: input_prefix !< Prefix of all input files, e.g. 'laf' for COSMO-DE analyses 76 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses 68 77 CHARACTER(LEN=SNAME) :: soil_prefix !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses 69 CHARACTER(LEN=SNAME) :: radiation_prefix !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses70 78 CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses 71 79 80 CHARACTER(LEN=SNAME) :: averaging_mode 72 81 CHARACTER(LEN=SNAME) :: bc_mode 73 82 CHARACTER(LEN=SNAME) :: ic_mode … … 77 86 REAL(dp) :: ug 78 87 REAL(dp) :: vg 79 REAL(dp) :: z0 !< Elevation of the PALM-4U domain above sea level [m] 88 REAL(dp) :: z0 !< Elevation of the PALM-4U domain above sea level [m] 89 REAL(dp) :: averaging_angle !< latitudal and longitudal width of averaging regions [deg] 90 91 92 LOGICAL :: debug 93 LOGICAL :: p0_is_set 94 LOGICAL :: ug_is_set 95 LOGICAL :: vg_is_set 96 LOGICAL :: flow_prefix_is_set !< 97 LOGICAL :: input_prefix_is_set !< 98 LOGICAL :: radiation_prefix_is_set !< 99 LOGICAL :: soil_prefix_is_set !< 100 LOGICAL :: soilmoisture_prefix_is_set !< 80 101 END TYPE inifor_config 81 102 82 103 TYPE grid_definition 83 104 CHARACTER(LEN=SNAME) :: name(3) !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/) 105 CHARACTER(LEN=SNAME) :: kind !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/) 106 INTEGER :: k_min !< Index of lowest PALM grid level that is not cut by local COSMO orography; vertically separates interpolation and extrapolation region. 84 107 INTEGER :: nx !< number of gridpoints in the first dimension 85 108 INTEGER :: ny !< number of gridpoints in the second dimension 86 INTEGER :: nz !< number of gridpoints in the third dimension 109 INTEGER :: nz !< number of gridpoints in the third dimension, used for PALM points 110 INTEGER :: nlev !< number of COSMO grid levels 111 INTEGER :: n_columns !< number of averaging columns of the source grid 87 112 INTEGER, ALLOCATABLE :: ii(:,:,:) !< Given a point (i,j,k) in the PALM-4U grid, ii(i,j,l) gives the x index of the l'th horizontl neighbour on the COSMO-DE grid. 88 113 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. 89 114 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. 115 INTEGER, ALLOCATABLE :: iii(:) !< profile averaging neighbour indices 116 INTEGER, ALLOCATABLE :: jjj(:) !< profile averaging neighbour indices 117 INTEGER, ALLOCATABLE :: kkk(:,:,:) !< indices of vertical interpolation neightbours, kkk(<source column>, <PALM k level>, <neighbour index>) 90 118 REAL(dp) :: lx !< domain length in the first dimension [m] 91 119 REAL(dp) :: ly !< domain length in the second dimension [m] … … 97 125 REAL(dp), POINTER :: z(:) !< coordinates of cell centers in z direction [m] 98 126 REAL(dp), ALLOCATABLE :: h(:,:,:) !< heights grid point for intermediate grids [m] 127 REAL(dp), POINTER :: cosmo_h(:,:,:)!< pointer to appropriate COSMO level heights (scalar/w) [m] 99 128 REAL(dp), POINTER :: hhl(:,:,:) !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from 100 129 REAL(dp), POINTER :: hfl(:,:,:) !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl … … 115 144 REAL(dp), ALLOCATABLE :: w_horiz(:,:,:) !< weights for bilinear horizontal interpolation 116 145 REAL(dp), ALLOCATABLE :: w_verti(:,:,:,:) !< weights for linear vertical interpolation 146 REAL(dp), ALLOCATABLE :: w(:,:,:) !< vertical interpolation weights, w(<source_column>, <PALM k level>, <neighbour index>) [-] 117 147 END TYPE grid_definition 118 148 … … 150 180 CHARACTER(LEN=SNAME) :: task !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile' 151 181 LOGICAL :: to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed 182 LOGICAL :: is_internal = .FALSE. !< INIFOR flag indicating whether variable shall be written to netCDF file (.FALSE.) or kept for later (.TRUE.) 152 183 LOGICAL :: is_read = .FALSE. !< INIFOR flag indicating whether variable has been read 153 LOGICAL :: is_upside_down = .FALSE. !< INIFOR flag indicating whether v ariable shall be processed184 LOGICAL :: is_upside_down = .FALSE. !< INIFOR flag indicating whether vertical dimension is reversed (typically the case with COSMO-DE atmospheric fields) 154 185 TYPE(grid_definition), POINTER :: grid !< Pointer to the corresponding output grid 155 186 TYPE(grid_definition), POINTER :: intermediate_grid !< Pointer to the corresponding intermediate grid 187 TYPE(grid_definition), POINTER :: averaging_grid !< Pointer to the corresponding intermediate grid 156 188 END TYPE nc_var 157 189 … … 159 191 TYPE io_group !< Input/Output group, groups together output variabels that share their input variables. For instance, all boundary surfaces and initialization fields of the potential temperature are base on T and p. 160 192 INTEGER :: nt !< maximum number of output time steps across all output variables 161 INTEGER :: nv !< number of output variables 193 INTEGER :: nv !< number of netCDF output variables 194 INTEGER :: n_inputs !< number of input variables 195 INTEGER :: n_output_quantities !< number of physical quantities required for computing netCDF output variables 162 196 CHARACTER(LEN=SNAME) :: kind !< kind of I/O group 163 197 CHARACTER(LEN=PATH), ALLOCATABLE :: in_files(:) !< list of nt input files -
palm/trunk/UTIL/inifor/src/util.f90
r3183 r3395 26 26 ! ----------------- 27 27 ! $Id$ 28 ! New routines for computing potential temperature and moist air density 29 ! Increased number of digits in real-to-str conversion 30 ! 31 ! 3183 2018-07-27 14:25:55Z suehring 28 32 ! Improved real-to-string conversion 29 33 ! … … 48 52 USE defs, & 49 53 ONLY : dp, PI, DATE, SNAME 54 USE types, & 55 ONLY : grid_definition 50 56 51 57 IMPLICIT NONE … … 281 287 282 288 289 !------------------------------------------------------------------------------! 290 ! Description: 291 ! ------------ 292 !> Converts the absolute temperature to the potential temperature in place using 293 !> the identity a^b = e^(b ln(a)). 294 !> 295 !> theta = T * (p_ref/p)^(R/c_p) = T * e^( R/c_p * ln(p_ref/p) ) 296 !------------------------------------------------------------------------------! 297 SUBROUTINE potential_temperature(t, p, p_ref, r, cp) 298 REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: t 299 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: p 300 REAL(dp), INTENT(IN) :: p_ref, r, cp 301 REAL(dp) :: rcp 302 303 rcp = r/cp 304 t(:,:,:) = t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) ) 305 306 END SUBROUTINE potential_temperature 307 308 309 !------------------------------------------------------------------------------! 310 ! Description: 311 ! ------------ 312 !> Compute the density in place of the given temperature (t_rho). 313 !------------------------------------------------------------------------------! 314 SUBROUTINE moist_density(t_rho, p, qv, rd, rv) 315 REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: t_rho 316 REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: p, qv 317 REAL(dp), INTENT(IN) :: rd, rv 318 319 t_rho(:,:,:) = p(:,:,:) / ( & 320 (rv * qv(:,:,:) + rd * (1.0_dp - qv(:,:,:))) * t_rho(:,:,:) & 321 ) 322 323 END SUBROUTINE moist_density 324 325 283 326 ! Convert a real number to a string in scientific notation 284 327 ! showing four significant digits. … … 298 341 299 342 300 CHARACTER(LEN=1 2) FUNCTION real_to_str_f(val)343 CHARACTER(LEN=16) FUNCTION real_to_str_f(val) 301 344 302 345 REAL(dp), INTENT(IN) :: val 303 346 304 WRITE(real_to_str_f, '(F1 2.4)') val347 WRITE(real_to_str_f, '(F16.8)') val 305 348 real_to_str_f = ADJUSTL(real_to_str_f) 306 349
Note: See TracChangeset
for help on using the changeset viewer.
