Changeset 3866 for palm/trunk/UTIL/inifor/src/inifor_grid.f90
- Timestamp:
- Apr 5, 2019 2:25:01 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/inifor/src/inifor_grid.f90
r3802 r3866 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Use PALM's working precision 29 ! Catch errors while reading namelists 30 ! Improved coding style 31 ! 32 ! 33 ! 3802 2019-03-17 13:33:42Z raasch 28 34 ! unused variable removed 29 35 ! … … 131 137 USE inifor_control 132 138 USE inifor_defs, & 133 ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,&139 ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, & 134 140 SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, & 135 141 BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM, & 136 RHO_L, OMEGA, HECTO 142 RHO_L, OMEGA, HECTO, wp, iwp 137 143 USE inifor_io, & 138 144 ONLY: get_cosmo_grid, get_netcdf_attribute, get_netcdf_dim_vector, & … … 140 146 get_input_file_list, validate_config 141 147 USE inifor_transform, & 142 ONLY: average_2d, rotate_to_cosmo, find_horizontal_neighbours, &148 ONLY: average_2d, rotate_to_cosmo, find_horizontal_neighbours, & 143 149 compute_horizontal_interp_weights, & 144 150 find_vertical_neighbours_and_weights_interp, & … … 156 162 SAVE 157 163 158 REAL( dp) :: averaging_angle = 0.0_dp !< latitudal and longitudal width of averaging regions [rad]159 REAL( dp) :: averaging_width_ns = 0.0_dp !< longitudal width of averaging regions [m]160 REAL( dp) :: averaging_width_ew = 0.0_dp !< latitudal width of averaging regions [m]161 REAL( dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad]162 REAL( dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad]163 REAL( dp) :: lambda_n = 0.0_dp !< longitude of rotaded pole of COSMO-DE grid [rad]164 REAL( dp) :: phi_c = 0.0_dp !< rotated-grid latitude of the center of the PALM domain [rad]165 REAL( dp) :: lambda_c = 0.0_dp !< rotated-grid longitude of the centre of the PALM domain [rad]166 REAL( dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad]167 REAL( dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad]168 REAL( dp) :: lam_centre = 0.0_dp !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]169 REAL( dp) :: phi_centre = 0.0_dp !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]170 REAL( dp) :: lam_east = 0.0_dp !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]171 REAL( dp) :: lam_west = 0.0_dp !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]172 REAL( dp) :: phi_north = 0.0_dp !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]173 REAL( dp) :: phi_south = 0.0_dp !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]174 REAL( dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug175 REAL( dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m]176 REAL( dp) :: dy = 0.0_dp !< PALM-4U grid spacing in y direction [m]177 REAL( dp) :: dz(10) = -1.0_dp !< PALM-4U grid spacing in z direction [m]178 REAL( dp) :: dz_max = 1000.0_dp !< maximum vertical grid spacing [m]179 REAL( dp) :: dz_stretch_factor = 1.08_dp !< factor for vertical grid stretching [m]180 REAL( dp) :: dz_stretch_level = -9999999.9_dp!< height above which the vertical grid will be stretched [m]181 REAL( dp) :: dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter182 REAL( dp) :: dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter183 REAL( dp) :: dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter184 REAL( dp) :: dxi = 0.0_dp !< inverse PALM-4U grid spacing in x direction [m^-1]185 REAL( dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1]186 REAL( dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1]187 REAL( dp) :: f3 = 0.0_dp !< Coriolis parameter188 REAL( dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m]189 REAL( dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m]190 REAL( dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa]191 REAL( dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m]192 REAL( dp) :: y0 = 0.0_dp !< y coordinate of PALM-4U Earth tangent [m]193 REAL( dp) :: z0 = 0.0_dp !< Elevation of the PALM-4U domain above sea level [m]194 REAL( dp) :: z_top = 0.0_dp !< height of the scalar top boundary [m]195 REAL( dp) :: zw_top = 0.0_dp !< height of the vertical velocity top boundary [m]196 REAL( dp) :: lonmin_cosmo = 0.0_dp !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]197 REAL( dp) :: lonmax_cosmo = 0.0_dp !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]198 REAL( dp) :: latmin_cosmo = 0.0_dp !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]199 REAL( dp) :: latmax_cosmo = 0.0_dp !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]200 REAL( dp) :: lonmin_palm = 0.0_dp !< Minimunm longitude of PALM grid [COSMO rotated-pole rad]201 REAL( dp) :: lonmax_palm = 0.0_dp !< Maximum longitude of PALM grid [COSMO rotated-pole rad]202 REAL( dp) :: latmin_palm = 0.0_dp !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]203 REAL( dp) :: latmax_palm = 0.0_dp !< Maximum latitude of PALM grid [COSMO rotated-pole rad]204 REAL( dp) :: lonmin_tot = 0.0_dp !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]205 REAL( dp) :: lonmax_tot = 0.0_dp !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]206 REAL( dp) :: latmin_tot = 0.0_dp !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]207 REAL( dp) :: latmax_tot = 0.0_dp !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]208 REAL( dp) :: latitude = 0.0_dp !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]209 REAL( dp) :: longitude = 0.0_dp !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]210 REAL( dp) :: origin_lat = 0.0_dp !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]211 REAL( dp) :: origin_lon = 0.0_dp !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]212 REAL( dp) :: rotation_angle = 0.0_dp !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]213 REAL( dp) :: end_time = 0.0_dp !< PALM-4U simulation time [s]214 215 REAL( dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file216 REAL( dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl217 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: depths !< COSMO-DE's TERRA-ML soil layer depths218 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth !< COSMO-DE's TERRA-ML soil layer thicknesses219 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth_rho_inv !< inverted soil water mass220 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlon !< longitudes of COSMO-DE's rotated-pole grid221 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlat !< latitudes of COSMO-DE's rotated-pole grid222 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: time !< output times223 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: x !< base palm grid x coordinate vector pointed to by grid_definitions224 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: xu !< base palm grid xu coordinate vector pointed to by grid_definitions225 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: y !< base palm grid y coordinate vector pointed to by grid_definitions226 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: yv !< base palm grid yv coordinate vector pointed to by grid_definitions227 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: z_column !< base palm grid z coordinate vector including the top boundary coordinate (entire column)228 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw_column !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)229 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: z !< base palm grid z coordinate vector pointed to by grid_definitions230 REAL( dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw !< base palm grid zw coordinate vector pointed to by grid_definitions231 232 INTEGER( hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: soiltyp!< COSMO-DE soil type map164 REAL(wp) :: averaging_angle = 0.0_wp !< latitudal and longitudal width of averaging regions [rad] 165 REAL(wp) :: averaging_width_ns = 0.0_wp !< longitudal width of averaging regions [m] 166 REAL(wp) :: averaging_width_ew = 0.0_wp !< latitudal width of averaging regions [m] 167 REAL(wp) :: phi_equat = 0.0_wp !< latitude of rotated equator of COSMO-DE grid [rad] 168 REAL(wp) :: phi_n = 0.0_wp !< latitude of rotated pole of COSMO-DE grid [rad] 169 REAL(wp) :: lambda_n = 0.0_wp !< longitude of rotaded pole of COSMO-DE grid [rad] 170 REAL(wp) :: phi_c = 0.0_wp !< rotated-grid latitude of the center of the PALM domain [rad] 171 REAL(wp) :: lambda_c = 0.0_wp !< rotated-grid longitude of the centre of the PALM domain [rad] 172 REAL(wp) :: phi_cn = 0.0_wp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] 173 REAL(wp) :: lambda_cn = 0.0_wp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] 174 REAL(wp) :: lam_centre = 0.0_wp !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 175 REAL(wp) :: phi_centre = 0.0_wp !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad] 176 REAL(wp) :: lam_east = 0.0_wp !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 177 REAL(wp) :: lam_west = 0.0_wp !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 178 REAL(wp) :: phi_north = 0.0_wp !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 179 REAL(wp) :: phi_south = 0.0_wp !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad] 180 REAL(wp) :: gam = 0.0_wp !< angle for working around phirot2phi/rlarot2rla bug 181 REAL(wp) :: dx = 0.0_wp !< PALM-4U grid spacing in x direction [m] 182 REAL(wp) :: dy = 0.0_wp !< PALM-4U grid spacing in y direction [m] 183 REAL(wp) :: dz(10) = -1.0_wp !< PALM-4U grid spacing in z direction [m] 184 REAL(wp) :: dz_max = 1000.0_wp !< maximum vertical grid spacing [m] 185 REAL(wp) :: dz_stretch_factor = 1.08_wp !< factor for vertical grid stretching [m] 186 REAL(wp) :: dz_stretch_level = -9999999.9_wp!< height above which the vertical grid will be stretched [m] 187 REAL(wp) :: dz_stretch_level_start(9) = -9999999.9_wp !< namelist parameter 188 REAL(wp) :: dz_stretch_level_end(9) = 9999999.9_wp !< namelist parameter 189 REAL(wp) :: dz_stretch_factor_array(9) = 1.08_wp !< namelist parameter 190 REAL(wp) :: dxi = 0.0_wp !< inverse PALM-4U grid spacing in x direction [m^-1] 191 REAL(wp) :: dyi = 0.0_wp !< inverse PALM-4U grid spacing in y direction [m^-1] 192 REAL(wp) :: dzi = 0.0_wp !< inverse PALM-4U grid spacing in z direction [m^-1] 193 REAL(wp) :: f3 = 0.0_wp !< Coriolis parameter 194 REAL(wp) :: lx = 0.0_wp !< PALM-4U domain size in x direction [m] 195 REAL(wp) :: ly = 0.0_wp !< PALM-4U domain size in y direction [m] 196 REAL(wp) :: p0 = 0.0_wp !< PALM-4U surface pressure, at z0 [Pa] 197 REAL(wp) :: x0 = 0.0_wp !< x coordinate of PALM-4U Earth tangent [m] 198 REAL(wp) :: y0 = 0.0_wp !< y coordinate of PALM-4U Earth tangent [m] 199 REAL(wp) :: z0 = 0.0_wp !< Elevation of the PALM-4U domain above sea level [m] 200 REAL(wp) :: z_top = 0.0_wp !< height of the scalar top boundary [m] 201 REAL(wp) :: zw_top = 0.0_wp !< height of the vertical velocity top boundary [m] 202 REAL(wp) :: lonmin_cosmo = 0.0_wp !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad] 203 REAL(wp) :: lonmax_cosmo = 0.0_wp !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad] 204 REAL(wp) :: latmin_cosmo = 0.0_wp !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad] 205 REAL(wp) :: latmax_cosmo = 0.0_wp !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad] 206 REAL(wp) :: lonmin_palm = 0.0_wp !< Minimunm longitude of PALM grid [COSMO rotated-pole rad] 207 REAL(wp) :: lonmax_palm = 0.0_wp !< Maximum longitude of PALM grid [COSMO rotated-pole rad] 208 REAL(wp) :: latmin_palm = 0.0_wp !< Minimunm latitude of PALM grid [COSMO rotated-pole rad] 209 REAL(wp) :: latmax_palm = 0.0_wp !< Maximum latitude of PALM grid [COSMO rotated-pole rad] 210 REAL(wp) :: lonmin_tot = 0.0_wp !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad] 211 REAL(wp) :: lonmax_tot = 0.0_wp !< Maximum longitude of required COSMO data [COSMO rotated-pole rad] 212 REAL(wp) :: latmin_tot = 0.0_wp !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad] 213 REAL(wp) :: latmax_tot = 0.0_wp !< Maximum latitude of required COSMO data [COSMO rotated-pole rad] 214 REAL(wp) :: latitude = 0.0_wp !< geographical latitude of the PALM-4U origin, from inipar namelist [deg] 215 REAL(wp) :: longitude = 0.0_wp !< geographical longitude of the PALM-4U origin, from inipar namelist [deg] 216 REAL(wp) :: origin_lat = 0.0_wp !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg] 217 REAL(wp) :: origin_lon = 0.0_wp !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg] 218 REAL(wp) :: rotation_angle = 0.0_wp !< clockwise angle the PALM-4U north is rotated away from geographical north [deg] 219 REAL(wp) :: end_time = 0.0_wp !< PALM-4U simulation time [s] 220 221 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file 222 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl 223 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: depths !< COSMO-DE's TERRA-ML soil layer depths 224 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth !< COSMO-DE's TERRA-ML soil layer thicknesses 225 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth_rho_inv !< inverted soil water mass 226 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: rlon !< longitudes of COSMO-DE's rotated-pole grid 227 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: rlat !< latitudes of COSMO-DE's rotated-pole grid 228 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: time !< output times 229 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: x !< base palm grid x coordinate vector pointed to by grid_definitions 230 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: xu !< base palm grid xu coordinate vector pointed to by grid_definitions 231 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: y !< base palm grid y coordinate vector pointed to by grid_definitions 232 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: yv !< base palm grid yv coordinate vector pointed to by grid_definitions 233 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: z_column !< base palm grid z coordinate vector including the top boundary coordinate (entire column) 234 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: zw_column !< base palm grid zw coordinate vector including the top boundary coordinate (entire column) 235 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: z !< base palm grid z coordinate vector pointed to by grid_definitions 236 REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: zw !< base palm grid zw coordinate vector pointed to by grid_definitions 237 238 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: soiltyp !< COSMO-DE soil type map 233 239 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 234 240 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 241 INTEGER :: iostat !< return status of READ statement 235 242 INTEGER :: nt !< number of output time steps 236 243 INTEGER :: nx !< number of PALM-4U grid points in x direction … … 355 362 !> interplation grids later in setup_grids(). 356 363 !------------------------------------------------------------------------------! 357 364 SUBROUTINE setup_parameters() 358 365 359 366 ! … … 361 368 ! Section 1: Define default parameters 362 369 !------------------------------------------------------------------------------ 363 cfg %start_date = '2013072100'364 365 366 367 368 ! 369 !-- 370 origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation371 origin_lon = 13.082744_dp * TO_RADIANS372 cfg % z0 = 35.0_dp373 374 ! 375 !-- 376 cfg % ug = 0.0_dp377 cfg % vg = 0.0_dp378 cfg %p0 = P_SL379 380 ! 381 !-- 382 383 384 385 386 387 388 389 390 cfg %flow_prefix = input_prefix391 cfg %input_prefix = input_prefix392 cfg %soil_prefix = input_prefix393 cfg %radiation_prefix = input_prefix394 cfg %soilmoisture_prefix = input_prefix395 396 397 398 399 400 401 cfg %debug = .FALSE.402 cfg % averaging_angle = 2.0_dp370 cfg%start_date = '2013072100' 371 end_hour = 2 372 start_hour_soil = -2 373 start_hour_soilmoisture = - (4 * 7 * 24) - 2 374 375 ! 376 !-- Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g) 377 origin_lat = 52.325079_wp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation 378 origin_lon = 13.082744_wp * TO_RADIANS 379 cfg%z0 = 35.0_wp 380 381 ! 382 !-- Default atmospheric parameters 383 cfg%ug = 0.0_wp 384 cfg%vg = 0.0_wp 385 cfg%p0 = P_SL 386 387 ! 388 !-- Parameters for file names 389 start_hour_flow = 0 390 start_hour_soil = 0 391 start_hour_radiation = 0 392 start_hour_soilmoisture = start_hour_flow - 2 393 end_hour_soilmoisture = start_hour_flow 394 step_hour = FORCING_STEP 395 396 input_prefix = 'laf' 397 cfg%flow_prefix = input_prefix 398 cfg%input_prefix = input_prefix 399 cfg%soil_prefix = input_prefix 400 cfg%radiation_prefix = input_prefix 401 cfg%soilmoisture_prefix = input_prefix 402 403 flow_suffix = '-flow' 404 soil_suffix = '-soil' 405 radiation_suffix = '-rad' 406 soilmoisture_suffix = '-soilmoisture' 407 408 cfg%debug = .FALSE. 409 cfg%averaging_angle = 2.0_wp 403 410 ! 404 411 !------------------------------------------------------------------------------ … … 407 414 408 415 ! 409 !-- Set default paths and modes 410 cfg % input_path = './' 411 cfg % hhl_file = '' 412 cfg % soiltyp_file = '' 413 cfg % namelist_file = './namelist' 414 cfg % static_driver_file = '' 415 cfg % output_file = './palm-4u-input.nc' 416 cfg % ic_mode = 'volume' 417 cfg % bc_mode = 'real' 418 cfg % averaging_mode = 'level' 419 420 ! 421 !-- Overwrite defaults with user configuration 422 CALL parse_command_line_arguments( cfg ) 423 CALL report('main_loop', 'Running INIFOR version ' // VERSION) 424 425 flow_prefix = TRIM(cfg % input_prefix) 426 radiation_prefix = TRIM(cfg % input_prefix) 427 soil_prefix = TRIM(cfg % input_prefix) 428 soilmoisture_prefix = TRIM(cfg % input_prefix) 429 IF (cfg % flow_prefix_is_set) flow_prefix = TRIM(cfg % flow_prefix) 430 IF (cfg % radiation_prefix_is_set) radiation_prefix = TRIM(cfg % radiation_prefix) 431 IF (cfg % soil_prefix_is_set) soil_prefix = TRIM(cfg % soil_prefix) 432 IF (cfg % soilmoisture_prefix_is_set) soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix) 433 434 output_file % name = cfg % output_file 435 z0 = cfg % z0 436 p0 = cfg % p0 437 438 init_variables_required = .TRUE. 439 boundary_variables_required = TRIM( cfg % bc_mode ) == 'real' 440 ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal' 441 surface_forcing_required = .FALSE. 442 443 IF ( ls_forcing_variables_required ) THEN 444 message = "Averaging of large-scale forcing profiles " // & 445 "has not been implemented, yet." 446 CALL inifor_abort('setup_parameters', message) 447 ENDIF 448 449 ! 450 !-- Set default file paths, if not specified by user. 451 CALL normalize_path(cfg % input_path) 452 IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' 453 IF (TRIM(cfg % soiltyp_file) == '') cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc' 454 455 CALL validate_config( cfg ) 456 457 CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode)) 458 CALL report('setup_parameters', " forcing mode: " // TRIM(cfg % bc_mode)) 459 CALL report('setup_parameters', " averaging mode: " // TRIM(cfg % averaging_mode)) 460 CALL report('setup_parameters', " averaging angle: " // real_to_str(cfg % averaging_angle)) 461 CALL report('setup_parameters', " data path: " // TRIM(cfg % input_path)) 462 CALL report('setup_parameters', " hhl file: " // TRIM(cfg % hhl_file)) 463 CALL report('setup_parameters', " soiltyp file: " // TRIM(cfg % soiltyp_file)) 464 CALL report('setup_parameters', " namelist file: " // TRIM(cfg % namelist_file)) 465 CALL report('setup_parameters', " output data file: " // TRIM(output_file % name)) 466 IF (cfg % debug ) CALL report('setup_parameters', " debugging mode: enabled") 467 468 CALL run_control('time', 'init') 469 ! 470 !-- Read in namelist parameters 471 OPEN(10, FILE=cfg % namelist_file) 472 READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz 473 READ(10, NML=d3par) ! end_time 416 !-- Set default paths and modes 417 cfg%input_path = './' 418 cfg%hhl_file = '' 419 cfg%soiltyp_file = '' 420 cfg%namelist_file = './namelist' 421 cfg%static_driver_file = '' 422 cfg%output_file = './palm-4u-input.nc' 423 cfg%ic_mode = 'volume' 424 cfg%bc_mode = 'real' 425 cfg%averaging_mode = 'level' 426 427 ! 428 !-- Overwrite defaults with user configuration 429 CALL parse_command_line_arguments( cfg ) 430 CALL report('main_loop', 'Running INIFOR version ' // VERSION) 431 432 flow_prefix = TRIM(cfg%input_prefix) 433 radiation_prefix = TRIM(cfg%input_prefix) 434 soil_prefix = TRIM(cfg%input_prefix) 435 soilmoisture_prefix = TRIM(cfg%input_prefix) 436 IF (cfg%flow_prefix_is_set) flow_prefix = TRIM(cfg%flow_prefix) 437 IF (cfg%radiation_prefix_is_set) radiation_prefix = TRIM(cfg%radiation_prefix) 438 IF (cfg%soil_prefix_is_set) soil_prefix = TRIM(cfg%soil_prefix) 439 IF (cfg%soilmoisture_prefix_is_set) soilmoisture_prefix = TRIM(cfg%soilmoisture_prefix) 440 441 output_file%name = cfg%output_file 442 z0 = cfg%z0 443 p0 = cfg%p0 444 445 init_variables_required = .TRUE. 446 boundary_variables_required = TRIM( cfg%bc_mode ) == 'real' 447 ls_forcing_variables_required = TRIM( cfg%bc_mode ) == 'ideal' 448 surface_forcing_required = .FALSE. 449 450 IF ( ls_forcing_variables_required ) THEN 451 message = "Averaging of large-scale forcing profiles " // & 452 "has not been implemented, yet." 453 CALL inifor_abort('setup_parameters', message) 454 ENDIF 455 456 ! 457 !-- Set default file paths, if not specified by user. 458 CALL normalize_path(cfg%input_path) 459 IF (TRIM(cfg%hhl_file) == '') cfg%hhl_file = TRIM(cfg%input_path) // 'hhl.nc' 460 IF (TRIM(cfg%soiltyp_file) == '') cfg%soiltyp_file = TRIM(cfg%input_path) // 'soil.nc' 461 462 CALL validate_config( cfg ) 463 464 CALL report('setup_parameters', "initialization mode: " // TRIM(cfg%ic_mode)) 465 CALL report('setup_parameters', " forcing mode: " // TRIM(cfg%bc_mode)) 466 CALL report('setup_parameters', " averaging mode: " // TRIM(cfg%averaging_mode)) 467 CALL report('setup_parameters', " averaging angle: " // real_to_str(cfg%averaging_angle)) 468 CALL report('setup_parameters', " data path: " // TRIM(cfg%input_path)) 469 CALL report('setup_parameters', " hhl file: " // TRIM(cfg%hhl_file)) 470 CALL report('setup_parameters', " soiltyp file: " // TRIM(cfg%soiltyp_file)) 471 CALL report('setup_parameters', " namelist file: " // TRIM(cfg%namelist_file)) 472 CALL report('setup_parameters', " output data file: " // TRIM(output_file%name)) 473 IF (cfg%debug ) CALL report('setup_parameters', " debugging mode: enabled") 474 475 CALL log_runtime('time', 'init') 476 ! 477 !-- Read in namelist parameters 478 OPEN(10, FILE=cfg%namelist_file, STATUS='old') 479 READ(10, NML=inipar, IOSTAT=iostat) ! nx, ny, nz, dx, dy, dz 480 IF ( iostat > 0 ) THEN 481 message = "Failed to read namelist 'inipar' from file '" // & 482 TRIM( cfg%namelist_file ) // "'. " 483 CALL inifor_abort( 'setup_parameters', message ) 474 484 CLOSE(10) 475 CALL run_control('time', 'read') 476 477 end_hour = CEILING( end_time / 3600.0 * step_hour ) 478 479 ! 480 !-- Generate input file lists 481 CALL get_input_file_list( & 482 cfg % start_date, start_hour_flow, end_hour, step_hour, & 483 cfg % input_path, flow_prefix, flow_suffix, flow_files) 484 CALL get_input_file_list( & 485 cfg % start_date, start_hour_soil, end_hour, step_hour, & 486 cfg % input_path, soil_prefix, soil_suffix, soil_files) 487 CALL get_input_file_list( & 488 cfg % start_date, start_hour_radiation, end_hour, step_hour, & 489 cfg % input_path, radiation_prefix, radiation_suffix, radiation_files, nocheck=.TRUE.) 490 CALL get_input_file_list( & 491 cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 492 cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files, nocheck=.TRUE.) 485 ENDIF 486 487 READ(10, NML=d3par, IOSTAT=iostat) ! end_time 488 IF ( iostat > 0 ) THEN 489 message = "Failed to read namelist 'd3par' from file '" // & 490 TRIM( cfg%namelist_file ) // "'. " 491 CALL inifor_abort( 'setup_parameters', message ) 492 CLOSE(10) 493 ENDIF 494 CLOSE(10) 495 496 CALL log_runtime('time', 'read') 497 498 end_hour = CEILING( end_time / 3600.0 * step_hour ) 499 500 ! 501 !-- Generate input file lists 502 CALL get_input_file_list( & 503 cfg%start_date, start_hour_flow, end_hour, step_hour, & 504 cfg%input_path, flow_prefix, flow_suffix, flow_files) 505 CALL get_input_file_list( & 506 cfg%start_date, start_hour_soil, end_hour, step_hour, & 507 cfg%input_path, soil_prefix, soil_suffix, soil_files) 508 CALL get_input_file_list( & 509 cfg%start_date, start_hour_radiation, end_hour, step_hour, & 510 cfg%input_path, radiation_prefix, radiation_suffix, radiation_files, nocheck=.TRUE.) 511 CALL get_input_file_list( & 512 cfg%start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & 513 cfg%input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files, nocheck=.TRUE.) 493 514 494 515 ! … … 506 527 507 528 508 CALL run_control('time', 'init')509 ! 510 !-- 511 cosmo_var %name = 'SOILTYP'512 CALL get_netcdf_variable(cfg %soiltyp_file, cosmo_var, soiltyp)513 514 515 IF (TRIM(cfg %static_driver_file) .NE. '') THEN516 517 origin_lon = get_netcdf_attribute(cfg %static_driver_file, 'origin_lon')518 origin_lat = get_netcdf_attribute(cfg %static_driver_file, 'origin_lat')519 520 521 // TRIM(cfg %static_driver_file) // "'"522 523 524 525 526 527 528 529 530 // TRIM(cfg %namelist_file) // "'"531 532 533 534 535 536 537 538 CALL run_control('time', 'read')539 540 541 542 543 544 545 529 CALL log_runtime('time', 'init') 530 ! 531 !-- Read COSMO soil type map 532 cosmo_var%name = 'SOILTYP' 533 CALL get_netcdf_variable(cfg%soiltyp_file, cosmo_var, soiltyp) 534 535 message = 'Reading PALM-4U origin from' 536 IF (TRIM(cfg%static_driver_file) .NE. '') THEN 537 538 origin_lon = get_netcdf_attribute(cfg%static_driver_file, 'origin_lon') 539 origin_lat = get_netcdf_attribute(cfg%static_driver_file, 'origin_lat') 540 541 message = TRIM(message) // " static driver file '" & 542 // TRIM(cfg%static_driver_file) // "'" 543 544 545 ELSE 546 547 origin_lon = longitude 548 origin_lat = latitude 549 550 message = TRIM(message) // " namlist file '" & 551 // TRIM(cfg%namelist_file) // "'" 552 553 ENDIF 554 origin_lon = origin_lon * TO_RADIANS 555 origin_lat = origin_lat * TO_RADIANS 556 557 CALL report('setup_parameters', message) 558 559 CALL log_runtime('time', 'read') 560 561 CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths, & 562 d_depth, d_depth_rho_inv, phi_n, lambda_n, & 563 phi_equat, & 564 lonmin_cosmo, lonmax_cosmo, & 565 latmin_cosmo, latmax_cosmo, & 566 nlon, nlat, nlev, ndepths ) 546 567 547 568 … … 550 571 !------------------------------------------------------------------------------ 551 572 ! 552 !-- 553 554 555 556 ! 557 !-- 558 x0 = 0.0_dp559 y0 = 0.0_dp560 561 ! 562 !-- 563 nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1564 565 CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time)566 output_file %time => time567 CALL run_control('time', 'init')568 569 ! 570 !-- 571 572 573 574 575 576 577 0.0_dp )578 579 ! 580 !-- 581 !-- 582 !-- 583 !-- 584 585 586 ! 587 !-- 588 !-- 589 !-- 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 CALL run_control('time', 'comp')573 !-- PALM-4U domain extents 574 lx = (nx+1) * dx 575 ly = (ny+1) * dy 576 577 ! 578 !-- PALM-4U point of Earth tangency 579 x0 = 0.0_wp 580 y0 = 0.0_wp 581 582 ! 583 !-- time vector 584 nt = CEILING(end_time / (step_hour * 3600.0_wp)) + 1 585 ALLOCATE( time(nt) ) 586 CALL linspace(0.0_wp, 3600.0_wp * (nt-1), time) 587 output_file%time => time 588 CALL log_runtime('time', 'init') 589 590 ! 591 !-- Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid 592 phi_c = TO_RADIANS * & 593 phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,& 594 phi_n * TO_DEGREES, lambda_n * TO_DEGREES ) 595 lambda_c = TO_RADIANS * & 596 rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,& 597 phi_n * TO_DEGREES, lambda_n * TO_DEGREES, & 598 0.0_wp ) 599 600 ! 601 !-- Set gamma according to whether PALM domain is in the northern or southern 602 !-- hemisphere of the COSMO rotated-pole system. Gamma assumes either the 603 !-- value 0 or PI and is needed to work around around a bug in the 604 !-- rotated-pole coordinate transformations. 605 gam = gamma_from_hemisphere(origin_lat, phi_equat) 606 607 ! 608 !-- Compute the north pole of the rotated-pole grid centred at the PALM-4U 609 !-- domain centre. The resulting (phi_cn, lambda_cn) are coordinates in 610 !-- COSMO-DE's rotated-pole grid. 611 phi_cn = phic_to_phin(phi_c) 612 lambda_cn = lamc_to_lamn(phi_c, lambda_c) 613 614 message = "PALM-4U origin:" // NEW_LINE('') // & 615 " lon (lambda) = " // & 616 TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //& 617 " lat (phi ) = " // & 618 TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //& 619 " lon (lambda) = " // & 620 TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & 621 " lat (phi ) = " // & 622 TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" 623 CALL report ('setup_parameters', message) 624 625 message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // & 626 " lon (lambda) = " // & 627 TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //& 628 " lat (phi ) = " // & 629 TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)" 630 CALL report ('setup_parameters', message) 631 632 message = "North pole of the rotated palm system:" // NEW_LINE(' ') // & 633 " lon (lambda) = " // & 634 TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & 635 " lat (phi ) = " // & 636 TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" 637 CALL report ('setup_parameters', message) 638 639 CALL log_runtime('time', 'comp') 619 640 620 641 !------------------------------------------------------------------------------ … … 625 646 !-- Compute coordiantes of the PALM centre in the source (COSMO) system 626 647 phi_centre = phirot2phi( & 627 phirot = project(0.5_ dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, &628 rlarot = project(0.5_ dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, &648 phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 649 rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 629 650 polphi = phi_cn * TO_DEGREES, & 630 651 polgam = gam * TO_DEGREES & … … 632 653 633 654 lam_centre = rlarot2rla( & 634 phirot = project(0.5_ dp*ly, y0, EARTH_RADIUS) * TO_DEGREES, &635 rlarot = project(0.5_ dp*lx, x0, EARTH_RADIUS) * TO_DEGREES, &655 phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES, & 656 rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES, & 636 657 polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, & 637 658 polgam = gam * TO_DEGREES & … … 647 668 ! 648 669 !-- Compute boundaries of the central averaging box 649 averaging_angle = cfg %averaging_angle * TO_RADIANS650 lam_east = lam_centre + 0.5_ dp * averaging_angle651 lam_west = lam_centre - 0.5_ dp * averaging_angle652 phi_north = phi_centre + 0.5_ dp * averaging_angle653 phi_south = phi_centre - 0.5_ dp * averaging_angle670 averaging_angle = cfg%averaging_angle * TO_RADIANS 671 lam_east = lam_centre + 0.5_wp * averaging_angle 672 lam_west = lam_centre - 0.5_wp * averaging_angle 673 phi_north = phi_centre + 0.5_wp * averaging_angle 674 phi_south = phi_centre - 0.5_wp * averaging_angle 654 675 averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS 655 676 averaging_width_ns = averaging_angle * EARTH_RADIUS … … 674 695 ! 675 696 !-- Coriolis parameter 676 f3 = 2.0_ dp * OMEGA * SIN( &697 f3 = 2.0_wp * OMEGA * SIN( & 677 698 TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,& 678 699 phi_n * TO_DEGREES, & … … 680 701 ) 681 702 682 703 END SUBROUTINE setup_parameters 683 704 684 705 … … 689 710 !> coordinates and interpolation weights 690 711 !------------------------------------------------------------------------------! 691 692 712 SUBROUTINE setup_grids() 713 CHARACTER :: interp_mode 693 714 694 715 !------------------------------------------------------------------------------ … … 696 717 !------------------------------------------------------------------------------ 697 718 ! 698 !-- 699 !-- 700 !-- 701 702 CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)703 CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y)704 705 706 707 708 709 710 711 712 713 ! 714 !-- 715 !-- 716 717 718 719 720 721 719 !-- palm x y z, we allocate the column to nz+1 in order to include the top 720 !-- scalar boundary. The interpolation grids will be associated with 721 !-- a shorter column that omits the top element. 722 ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) ) 723 CALL linspace(0.5_wp * dx, lx - 0.5_wp * dx, x) 724 CALL linspace(0.5_wp * dy, ly - 0.5_wp * dy, y) 725 CALL stretched_z(z_column, dz, dz_max=dz_max, & 726 dz_stretch_factor=dz_stretch_factor, & 727 dz_stretch_level=dz_stretch_level, & 728 dz_stretch_level_start=dz_stretch_level_start, & 729 dz_stretch_level_end=dz_stretch_level_end, & 730 dz_stretch_factor_array=dz_stretch_factor_array) 731 z(1:nz) = z_column(1:nz) 732 z_top = z_column(nz+1) 733 734 ! 735 !-- palm xu yv zw, compared to the scalar grid, velocity coordinates 736 !-- contain one element less. 737 ALLOCATE( xu(1:nx), yv(1:ny), zw(1:nz-1), zw_column(1:nz)) 738 CALL linspace(dx, lx - dx, xu) 739 CALL linspace(dy, ly - dy, yv) 740 CALL midpoints(z_column, zw_column) 741 zw(1:nz-1) = zw_column(1:nz-1) 742 zw_top = zw_column(nz) 722 743 723 744 … … 725 746 ! Section 1: Define initialization and boundary grids 726 747 !------------------------------------------------------------------------------ 727 CALL init_grid_definition('palm', grid=palm_grid, & 728 xmin=0.0_dp, xmax=lx, & 729 ymin=0.0_dp, ymax=ly, & 730 x0=x0, y0=y0, z0=z0, & 731 nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode) 732 733 ! 734 !-- Subtracting 1 because arrays will be allocated with nlon + 1 elements. 735 CALL init_grid_definition('cosmo-de', grid=cosmo_grid, & 736 xmin=lonmin_cosmo, xmax=lonmax_cosmo, & 737 ymin=latmin_cosmo, ymax=latmax_cosmo, & 738 x0=x0, y0=y0, z0=0.0_dp, & 739 nx=nlon-1, ny=nlat-1, nz=nlev-1) 740 741 ! 742 !-- Define intermediate grid. This is the same as palm_grid except with a 743 !-- much coarser vertical grid. The vertical levels are interpolated in each 744 !-- PALM column from COSMO's secondary levels. The main levels are then 745 !-- computed as the averages of the bounding secondary levels. 746 CALL init_grid_definition('palm intermediate', grid=palm_intermediate, & 747 xmin=0.0_dp, xmax=lx, & 748 ymin=0.0_dp, ymax=ly, & 749 x0=x0, y0=y0, z0=z0, & 750 nx=nx, ny=ny, nz=nlev-2) 751 752 CALL init_grid_definition('boundary', grid=u_initial_grid, & 748 CALL init_grid_definition('palm', grid=palm_grid, & 749 xmin=0.0_wp, xmax=lx, & 750 ymin=0.0_wp, ymax=ly, & 751 x0=x0, y0=y0, z0=z0, & 752 nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg%ic_mode) 753 754 ! 755 !-- Subtracting 1 because arrays will be allocated with nlon + 1 elements. 756 CALL init_grid_definition('cosmo-de', grid=cosmo_grid, & 757 xmin=lonmin_cosmo, xmax=lonmax_cosmo, & 758 ymin=latmin_cosmo, ymax=latmax_cosmo, & 759 x0=x0, y0=y0, z0=0.0_wp, & 760 nx=nlon-1, ny=nlat-1, nz=nlev-1) 761 762 ! 763 !-- Define intermediate grid. This is the same as palm_grid except with a 764 !-- much coarser vertical grid. The vertical levels are interpolated in each 765 !-- PALM column from COSMO's secondary levels. The main levels are then 766 !-- computed as the averages of the bounding secondary levels. 767 CALL init_grid_definition('palm intermediate', grid=palm_intermediate, & 768 xmin=0.0_wp, xmax=lx, & 769 ymin=0.0_wp, ymax=ly, & 770 x0=x0, y0=y0, z0=z0, & 771 nx=nx, ny=ny, nz=nlev-2) 772 773 CALL init_grid_definition('boundary', grid=u_initial_grid, & 774 xmin = dx, xmax = lx - dx, & 775 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 776 x0=x0, y0=y0, z0 = z0, & 777 nx = nx-1, ny = ny, nz = nz, & 778 z=z, ic_mode=cfg%ic_mode) 779 780 CALL init_grid_definition('boundary', grid=v_initial_grid, & 781 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 782 ymin = dy, ymax = ly - dy, & 783 x0=x0, y0=y0, z0 = z0, & 784 nx = nx, ny = ny-1, nz = nz, & 785 z=z, ic_mode=cfg%ic_mode) 786 787 CALL init_grid_definition('boundary', grid=w_initial_grid, & 788 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 789 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 790 x0=x0, y0=y0, z0 = z0, & 791 nx = nx, ny = ny, nz = nz-1, & 792 z=zw, ic_mode=cfg%ic_mode) 793 794 CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate, & 795 xmin = dx, xmax = lx - dx, & 796 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 797 x0=x0, y0=y0, z0 = z0, & 798 nx = nx-1, ny = ny, nz = nlev - 2) 799 800 CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate, & 801 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 802 ymin = dy, ymax = ly - dy, & 803 x0=x0, y0=y0, z0 = z0, & 804 nx = nx, ny = ny-1, nz = nlev - 2) 805 806 CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate, & 807 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 808 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 809 x0=x0, y0=y0, z0 = z0, & 810 nx = nx, ny = ny, nz = nlev - 1) 811 812 IF (boundary_variables_required) THEN 813 ! 814 !------------------------------------------------------------------------------ 815 ! Section 2: Define PALM-4U boundary grids 816 !------------------------------------------------------------------------------ 817 CALL init_grid_definition('boundary', grid=scalars_east_grid, & 818 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 819 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 820 x0=x0, y0=y0, z0 = z0, & 821 nx = 0, ny = ny, nz = nz, z=z) 822 823 CALL init_grid_definition('boundary', grid=scalars_west_grid, & 824 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 825 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 826 x0=x0, y0=y0, z0 = z0, & 827 nx = 0, ny = ny, nz = nz, z=z) 828 829 CALL init_grid_definition('boundary', grid=scalars_north_grid, & 830 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 831 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 832 x0=x0, y0=y0, z0 = z0, & 833 nx = nx, ny = 0, nz = nz, z=z) 834 835 CALL init_grid_definition('boundary', grid=scalars_south_grid, & 836 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 837 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 838 x0=x0, y0=y0, z0 = z0, & 839 nx = nx, ny = 0, nz = nz, z=z) 840 841 CALL init_grid_definition('boundary', grid=scalars_top_grid, & 842 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 843 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 844 x0=x0, y0=y0, z0 = z0, & 845 nx = nx, ny = ny, nz = 1, z=(/z_top/)) 846 847 CALL init_grid_definition('boundary', grid=u_east_grid, & 848 xmin = lx, xmax = lx, & 849 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 850 x0=x0, y0=y0, z0 = z0, & 851 nx = 0, ny = ny, nz = nz, z=z) 852 853 CALL init_grid_definition('boundary', grid=u_west_grid, & 854 xmin = 0.0_wp, xmax = 0.0_wp, & 855 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 856 x0=x0, y0=y0, z0 = z0, & 857 nx = 0, ny = ny, nz = nz, z=z) 858 859 CALL init_grid_definition('boundary', grid=u_north_grid, & 753 860 xmin = dx, xmax = lx - dx, & 754 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,&861 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 755 862 x0=x0, y0=y0, z0 = z0, & 756 nx = nx-1, ny = ny, nz = nz, & 757 z=z, ic_mode=cfg % ic_mode) 758 759 CALL init_grid_definition('boundary', grid=v_initial_grid, & 760 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 863 nx = nx-1, ny = 0, nz = nz, z=z) 864 865 CALL init_grid_definition('boundary', grid=u_south_grid, & 866 xmin = dx, xmax = lx - dx, & 867 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 868 x0=x0, y0=y0, z0 = z0, & 869 nx = nx-1, ny = 0, nz = nz, z=z) 870 871 CALL init_grid_definition('boundary', grid=u_top_grid, & 872 xmin = dx, xmax = lx - dx, & 873 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 874 x0=x0, y0=y0, z0 = z0, & 875 nx = nx-1, ny = ny, nz = 1, z=(/z_top/)) 876 877 CALL init_grid_definition('boundary', grid=v_east_grid, & 878 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 761 879 ymin = dy, ymax = ly - dy, & 762 880 x0=x0, y0=y0, z0 = z0, & 763 nx = nx, ny = ny-1, nz = nz, & 764 z=z, ic_mode=cfg % ic_mode) 765 766 CALL init_grid_definition('boundary', grid=w_initial_grid, & 767 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 768 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 881 nx = 0, ny = ny-1, nz = nz, z=z) 882 883 CALL init_grid_definition('boundary', grid=v_west_grid, & 884 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 885 ymin = dy, ymax = ly - dy, & 769 886 x0=x0, y0=y0, z0 = z0, & 770 nx = nx, ny = ny, nz = nz-1, & 771 z=zw, ic_mode=cfg % ic_mode) 772 773 CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate, & 887 nx = 0, ny = ny-1, nz = nz, z=z) 888 889 CALL init_grid_definition('boundary', grid=v_north_grid, & 890 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 891 ymin = ly, ymax = ly, & 892 x0=x0, y0=y0, z0 = z0, & 893 nx = nx, ny = 0, nz = nz, z=z) 894 895 CALL init_grid_definition('boundary', grid=v_south_grid, & 896 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 897 ymin = 0.0_wp, ymax = 0.0_wp, & 898 x0=x0, y0=y0, z0 = z0, & 899 nx = nx, ny = 0, nz = nz, z=z) 900 901 CALL init_grid_definition('boundary', grid=v_top_grid, & 902 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 903 ymin = dy, ymax = ly - dy, & 904 x0=x0, y0=y0, z0 = z0, & 905 nx = nx, ny = ny-1, nz = 1, z=(/z_top/)) 906 907 CALL init_grid_definition('boundary', grid=w_east_grid, & 908 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 909 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 910 x0=x0, y0=y0, z0 = z0, & 911 nx = 0, ny = ny, nz = nz - 1, z=zw) 912 913 CALL init_grid_definition('boundary', grid=w_west_grid, & 914 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 915 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 916 x0=x0, y0=y0, z0 = z0, & 917 nx = 0, ny = ny, nz = nz - 1, z=zw) 918 919 CALL init_grid_definition('boundary', grid=w_north_grid, & 920 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 921 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 922 x0=x0, y0=y0, z0 = z0, & 923 nx = nx, ny = 0, nz = nz - 1, z=zw) 924 925 CALL init_grid_definition('boundary', grid=w_south_grid, & 926 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 927 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 928 x0=x0, y0=y0, z0 = z0, & 929 nx = nx, ny = 0, nz = nz - 1, z=zw) 930 931 CALL init_grid_definition('boundary', grid=w_top_grid, & 932 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 933 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 934 x0=x0, y0=y0, z0 = z0, & 935 nx = nx, ny = ny, nz = 1, z=(/zw_top/)) 936 937 CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate, & 938 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 939 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 940 x0=x0, y0=y0, z0 = z0, & 941 nx = 0, ny = ny, nz = nlev - 2) 942 943 CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate, & 944 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 945 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 946 x0=x0, y0=y0, z0 = z0, & 947 nx = 0, ny = ny, nz = nlev - 2) 948 949 CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate, & 950 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 951 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 952 x0=x0, y0=y0, z0 = z0, & 953 nx = nx, ny = 0, nz = nlev - 2) 954 955 CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate, & 956 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 957 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 958 x0=x0, y0=y0, z0 = z0, & 959 nx = nx, ny = 0, nz = nlev - 2) 960 961 CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate, & 962 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 963 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 964 x0=x0, y0=y0, z0 = z0, & 965 nx = nx, ny = ny, nz = nlev - 2) 966 967 CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate, & 968 xmin = lx, xmax = lx, & 969 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 970 x0=x0, y0=y0, z0 = z0, & 971 nx = 0, ny = ny, nz = nlev - 2) 972 973 CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate, & 974 xmin = 0.0_wp, xmax = 0.0_wp, & 975 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 976 x0=x0, y0=y0, z0 = z0, & 977 nx = 0, ny = ny, nz = nlev - 2) 978 979 CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate, & 774 980 xmin = dx, xmax = lx - dx, & 775 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 981 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 982 x0=x0, y0=y0, z0 = z0, & 983 nx = nx-1, ny = 0, nz = nlev - 2) 984 985 CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate, & 986 xmin = dx, xmax = lx - dx, & 987 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 988 x0=x0, y0=y0, z0 = z0, & 989 nx = nx-1, ny = 0, nz = nlev - 2) 990 991 CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate, & 992 xmin = dx, xmax = lx - dx, & 993 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 776 994 x0=x0, y0=y0, z0 = z0, & 777 995 nx = nx-1, ny = ny, nz = nlev - 2) 778 996 779 CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate, & 780 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 997 CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate, & 998 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 999 ymin = dy, ymax = ly - dy, & 1000 x0=x0, y0=y0, z0 = z0, & 1001 nx = 0, ny = ny-1, nz = nlev - 2) 1002 1003 CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate, & 1004 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 1005 ymin = dy, ymax = ly - dy, & 1006 x0=x0, y0=y0, z0 = z0, & 1007 nx = 0, ny = ny-1, nz = nlev - 2) 1008 1009 CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate, & 1010 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 1011 ymin = ly, ymax = ly, & 1012 x0=x0, y0=y0, z0 = z0, & 1013 nx = nx, ny = 0, nz = nlev - 2) 1014 1015 CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate, & 1016 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 1017 ymin = 0.0_wp, ymax = 0.0_wp, & 1018 x0=x0, y0=y0, z0 = z0, & 1019 nx = nx, ny = 0, nz = nlev - 2) 1020 1021 CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate, & 1022 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 781 1023 ymin = dy, ymax = ly - dy, & 782 1024 x0=x0, y0=y0, z0 = z0, & 783 1025 nx = nx, ny = ny-1, nz = nlev - 2) 784 1026 785 CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate, & 786 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 787 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 1027 CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate, & 1028 xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx, & 1029 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 1030 x0=x0, y0=y0, z0 = z0, & 1031 nx = 0, ny = ny, nz = nlev - 1) 1032 1033 CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & 1034 xmin = -0.5_wp * dx, xmax = -0.5_wp * dx, & 1035 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 1036 x0=x0, y0=y0, z0 = z0, & 1037 nx = 0, ny = ny, nz = nlev - 1) 1038 1039 CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & 1040 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 1041 ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy, & 1042 x0=x0, y0=y0, z0 = z0, & 1043 nx = nx, ny = 0, nz = nlev - 1) 1044 1045 CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & 1046 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 1047 ymin = -0.5_wp * dy, ymax = -0.5_wp * dy, & 1048 x0=x0, y0=y0, z0 = z0, & 1049 nx = nx, ny = 0, nz = nlev - 1) 1050 1051 CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & 1052 xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx, & 1053 ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy, & 788 1054 x0=x0, y0=y0, z0 = z0, & 789 1055 nx = nx, ny = ny, nz = nlev - 1) 790 791 IF (boundary_variables_required) THEN 792 ! 793 !------------------------------------------------------------------------------ 794 ! Section 2: Define PALM-4U boundary grids 795 !------------------------------------------------------------------------------ 796 CALL init_grid_definition('boundary', grid=scalars_east_grid, & 797 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 798 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 799 x0=x0, y0=y0, z0 = z0, & 800 nx = 0, ny = ny, nz = nz, z=z) 801 802 CALL init_grid_definition('boundary', grid=scalars_west_grid, & 803 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 804 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 805 x0=x0, y0=y0, z0 = z0, & 806 nx = 0, ny = ny, nz = nz, z=z) 807 808 CALL init_grid_definition('boundary', grid=scalars_north_grid, & 809 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 810 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 811 x0=x0, y0=y0, z0 = z0, & 812 nx = nx, ny = 0, nz = nz, z=z) 813 814 CALL init_grid_definition('boundary', grid=scalars_south_grid, & 815 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 816 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 817 x0=x0, y0=y0, z0 = z0, & 818 nx = nx, ny = 0, nz = nz, z=z) 819 820 CALL init_grid_definition('boundary', grid=scalars_top_grid, & 821 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 822 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 823 x0=x0, y0=y0, z0 = z0, & 824 nx = nx, ny = ny, nz = 1, z=(/z_top/)) 825 826 CALL init_grid_definition('boundary', grid=u_east_grid, & 827 xmin = lx, xmax = lx, & 828 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 829 x0=x0, y0=y0, z0 = z0, & 830 nx = 0, ny = ny, nz = nz, z=z) 831 832 CALL init_grid_definition('boundary', grid=u_west_grid, & 833 xmin = 0.0_dp, xmax = 0.0_dp, & 834 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 835 x0=x0, y0=y0, z0 = z0, & 836 nx = 0, ny = ny, nz = nz, z=z) 837 838 CALL init_grid_definition('boundary', grid=u_north_grid, & 839 xmin = dx, xmax = lx - dx, & 840 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 841 x0=x0, y0=y0, z0 = z0, & 842 nx = nx-1, ny = 0, nz = nz, z=z) 843 844 CALL init_grid_definition('boundary', grid=u_south_grid, & 845 xmin = dx, xmax = lx - dx, & 846 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 847 x0=x0, y0=y0, z0 = z0, & 848 nx = nx-1, ny = 0, nz = nz, z=z) 849 850 CALL init_grid_definition('boundary', grid=u_top_grid, & 851 xmin = dx, xmax = lx - dx, & 852 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 853 x0=x0, y0=y0, z0 = z0, & 854 nx = nx-1, ny = ny, nz = 1, z=(/z_top/)) 855 856 CALL init_grid_definition('boundary', grid=v_east_grid, & 857 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 858 ymin = dy, ymax = ly - dy, & 859 x0=x0, y0=y0, z0 = z0, & 860 nx = 0, ny = ny-1, nz = nz, z=z) 861 862 CALL init_grid_definition('boundary', grid=v_west_grid, & 863 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 864 ymin = dy, ymax = ly - dy, & 865 x0=x0, y0=y0, z0 = z0, & 866 nx = 0, ny = ny-1, nz = nz, z=z) 867 868 CALL init_grid_definition('boundary', grid=v_north_grid, & 869 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 870 ymin = ly, ymax = ly, & 871 x0=x0, y0=y0, z0 = z0, & 872 nx = nx, ny = 0, nz = nz, z=z) 873 874 CALL init_grid_definition('boundary', grid=v_south_grid, & 875 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 876 ymin = 0.0_dp, ymax = 0.0_dp, & 877 x0=x0, y0=y0, z0 = z0, & 878 nx = nx, ny = 0, nz = nz, z=z) 879 880 CALL init_grid_definition('boundary', grid=v_top_grid, & 881 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 882 ymin = dy, ymax = ly - dy, & 883 x0=x0, y0=y0, z0 = z0, & 884 nx = nx, ny = ny-1, nz = 1, z=(/z_top/)) 885 886 CALL init_grid_definition('boundary', grid=w_east_grid, & 887 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 888 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 889 x0=x0, y0=y0, z0 = z0, & 890 nx = 0, ny = ny, nz = nz - 1, z=zw) 891 892 CALL init_grid_definition('boundary', grid=w_west_grid, & 893 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 894 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 895 x0=x0, y0=y0, z0 = z0, & 896 nx = 0, ny = ny, nz = nz - 1, z=zw) 897 898 CALL init_grid_definition('boundary', grid=w_north_grid, & 899 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 900 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 901 x0=x0, y0=y0, z0 = z0, & 902 nx = nx, ny = 0, nz = nz - 1, z=zw) 903 904 CALL init_grid_definition('boundary', grid=w_south_grid, & 905 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 906 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 907 x0=x0, y0=y0, z0 = z0, & 908 nx = nx, ny = 0, nz = nz - 1, z=zw) 909 910 CALL init_grid_definition('boundary', grid=w_top_grid, & 911 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 912 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 913 x0=x0, y0=y0, z0 = z0, & 914 nx = nx, ny = ny, nz = 1, z=(/zw_top/)) 915 916 CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate, & 917 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 918 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 919 x0=x0, y0=y0, z0 = z0, & 920 nx = 0, ny = ny, nz = nlev - 2) 921 922 CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate, & 923 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 924 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 925 x0=x0, y0=y0, z0 = z0, & 926 nx = 0, ny = ny, nz = nlev - 2) 927 928 CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate, & 929 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 930 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 931 x0=x0, y0=y0, z0 = z0, & 932 nx = nx, ny = 0, nz = nlev - 2) 933 934 CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate, & 935 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 936 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 937 x0=x0, y0=y0, z0 = z0, & 938 nx = nx, ny = 0, nz = nlev - 2) 939 940 CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate, & 941 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 942 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 943 x0=x0, y0=y0, z0 = z0, & 944 nx = nx, ny = ny, nz = nlev - 2) 945 946 CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate, & 947 xmin = lx, xmax = lx, & 948 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 949 x0=x0, y0=y0, z0 = z0, & 950 nx = 0, ny = ny, nz = nlev - 2) 951 952 CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate, & 953 xmin = 0.0_dp, xmax = 0.0_dp, & 954 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 955 x0=x0, y0=y0, z0 = z0, & 956 nx = 0, ny = ny, nz = nlev - 2) 957 958 CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate, & 959 xmin = dx, xmax = lx - dx, & 960 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 961 x0=x0, y0=y0, z0 = z0, & 962 nx = nx-1, ny = 0, nz = nlev - 2) 963 964 CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate, & 965 xmin = dx, xmax = lx - dx, & 966 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 967 x0=x0, y0=y0, z0 = z0, & 968 nx = nx-1, ny = 0, nz = nlev - 2) 969 970 CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate, & 971 xmin = dx, xmax = lx - dx, & 972 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 973 x0=x0, y0=y0, z0 = z0, & 974 nx = nx-1, ny = ny, nz = nlev - 2) 975 976 CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate, & 977 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 978 ymin = dy, ymax = ly - dy, & 979 x0=x0, y0=y0, z0 = z0, & 980 nx = 0, ny = ny-1, nz = nlev - 2) 981 982 CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate, & 983 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 984 ymin = dy, ymax = ly - dy, & 985 x0=x0, y0=y0, z0 = z0, & 986 nx = 0, ny = ny-1, nz = nlev - 2) 987 988 CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate, & 989 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 990 ymin = ly, ymax = ly, & 991 x0=x0, y0=y0, z0 = z0, & 992 nx = nx, ny = 0, nz = nlev - 2) 993 994 CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate, & 995 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 996 ymin = 0.0_dp, ymax = 0.0_dp, & 997 x0=x0, y0=y0, z0 = z0, & 998 nx = nx, ny = 0, nz = nlev - 2) 999 1000 CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate, & 1001 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 1002 ymin = dy, ymax = ly - dy, & 1003 x0=x0, y0=y0, z0 = z0, & 1004 nx = nx, ny = ny-1, nz = nlev - 2) 1005 1006 CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate, & 1007 xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & 1008 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 1009 x0=x0, y0=y0, z0 = z0, & 1010 nx = 0, ny = ny, nz = nlev - 1) 1011 1012 CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & 1013 xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & 1014 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 1015 x0=x0, y0=y0, z0 = z0, & 1016 nx = 0, ny = ny, nz = nlev - 1) 1017 1018 CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & 1019 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 1020 ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & 1021 x0=x0, y0=y0, z0 = z0, & 1022 nx = nx, ny = 0, nz = nlev - 1) 1023 1024 CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & 1025 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 1026 ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & 1027 x0=x0, y0=y0, z0 = z0, & 1028 nx = nx, ny = 0, nz = nlev - 1) 1029 1030 CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & 1031 xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & 1032 ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & 1033 x0=x0, y0=y0, z0 = z0, & 1034 nx = nx, ny = ny, nz = nlev - 1) 1035 ENDIF 1056 ENDIF 1036 1057 1037 1058 ! … … 1040 1061 !------------------------------------------------------------------------------ 1041 1062 1042 lonmin_palm = MINVAL(palm_intermediate %clon)1043 lonmax_palm = MAXVAL(palm_intermediate %clon)1044 latmin_palm = MINVAL(palm_intermediate %clat)1045 latmax_palm = MAXVAL(palm_intermediate %clat)1046 1047 1048 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1049 1050 1051 1052 1053 1054 x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, &1055 1056 1057 1058 1059 1060 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1061 1062 1063 1064 1065 1066 x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, &1067 1068 1069 1070 1071 1072 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1073 1074 1075 1076 1077 1078 1079 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1080 1081 1082 1083 1084 1085 1086 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1087 1088 1089 1090 1091 1092 1093 x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &1094 1095 1096 1097 1063 lonmin_palm = MINVAL(palm_intermediate%clon) 1064 lonmax_palm = MAXVAL(palm_intermediate%clon) 1065 latmin_palm = MINVAL(palm_intermediate%clat) 1066 latmax_palm = MAXVAL(palm_intermediate%clat) 1067 1068 CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid, & 1069 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1070 lonmin = lonmin_palm, lonmax = lonmax_palm, & 1071 latmin = latmin_palm, latmax = latmax_palm, & 1072 kind='scalar', name='averaged initial scalar') 1073 1074 CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid, & 1075 x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0, & 1076 lonmin = lonmin_palm, lonmax = lonmax_palm, & 1077 latmin = latmin_palm, latmax = latmax_palm, & 1078 kind='w', name='averaged initial w') 1079 1080 CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid, & 1081 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1082 lonmin = lam_west, lonmax = lam_east, & 1083 latmin = phi_south, latmax = phi_north, & 1084 kind='scalar', name='centre geostrophic scalar') 1085 1086 CALL init_averaging_grid(averaged_w_profile, cosmo_grid, & 1087 x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0, & 1088 lonmin = lam_west, lonmax = lam_east, & 1089 latmin = phi_south, latmax = phi_north, & 1090 kind='w', name='centre geostrophic w') 1091 1092 CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid, & 1093 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1094 lonmin = lam_west, lonmax = lam_east, & 1095 latmin = phi_centre - averaging_angle, & 1096 latmax = phi_centre, & 1097 kind='scalar', name='south geostrophic scalar') 1098 1099 CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid, & 1100 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1101 lonmin = lam_west, lonmax = lam_east, & 1102 latmin = phi_centre, & 1103 latmax = phi_centre + averaging_angle, & 1104 kind='scalar', name='north geostrophic scalar') 1105 1106 CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid, & 1107 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1108 lonmin = lam_centre - averaging_angle, & 1109 lonmax = lam_centre, & 1110 latmin = phi_south, latmax = phi_north, & 1111 kind='scalar', name='west geostrophic scalar') 1112 1113 CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid, & 1114 x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0, & 1115 lonmin = lam_centre, & 1116 lonmax = lam_centre + averaging_angle, & 1117 latmin = phi_south, latmax = phi_north, & 1118 kind='scalar', name='east geostrophic scalar') 1098 1119 1099 1120 … … 1102 1123 ! Section 4: Precompute neighbours and weights for interpolation 1103 1124 !------------------------------------------------------------------------------ 1104 interp_mode = 's' 1105 CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode) 1106 IF (boundary_variables_required) THEN 1107 CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode) 1108 CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode) 1109 CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode) 1110 CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode) 1111 CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode) 1112 ENDIF 1113 1114 interp_mode = 'u' 1115 CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 1116 IF (boundary_variables_required) THEN 1117 CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode) 1118 CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode) 1119 CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode) 1120 CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode) 1121 CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode) 1122 ENDIF 1123 1124 interp_mode = 'v' 1125 CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 1126 IF (boundary_variables_required) THEN 1127 CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode) 1128 CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode) 1129 CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode) 1130 CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode) 1131 CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode) 1132 ENDIF 1133 1134 interp_mode = 'w' 1135 CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) 1136 IF (boundary_variables_required) THEN 1137 CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode) 1138 CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode) 1139 CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode) 1140 CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode) 1141 CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode) 1142 ENDIF 1143 1144 IF (TRIM(cfg % ic_mode) == 'profile') THEN 1145 !TODO: remove this conditional if not needed. 1146 ENDIF 1147 1148 1149 END SUBROUTINE setup_grids 1125 interp_mode = 's' 1126 CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg%ic_mode) 1127 IF (boundary_variables_required) THEN 1128 CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode) 1129 CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode) 1130 CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode) 1131 CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode) 1132 CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode) 1133 ENDIF 1134 1135 interp_mode = 'u' 1136 CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode) 1137 IF (boundary_variables_required) THEN 1138 CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode) 1139 CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode) 1140 CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode) 1141 CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode) 1142 CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode) 1143 ENDIF 1144 1145 interp_mode = 'v' 1146 CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode) 1147 IF (boundary_variables_required) THEN 1148 CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode) 1149 CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode) 1150 CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode) 1151 CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode) 1152 CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode) 1153 ENDIF 1154 1155 interp_mode = 'w' 1156 CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode) 1157 IF (boundary_variables_required) THEN 1158 CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode) 1159 CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode) 1160 CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode) 1161 CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode) 1162 CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode) 1163 ENDIF 1164 1165 IF (TRIM(cfg%ic_mode) == 'profile') THEN 1166 !TODO: remove this conditional if not needed. 1167 ENDIF 1168 1169 END SUBROUTINE setup_grids 1150 1170 1151 1171 … … 1156 1176 !> vertical interpolation. 1157 1177 !------------------------------------------------------------------------------! 1158 1159 1160 1161 1162 1163 1164 1165 REAL(dp), DIMENSION(:), POINTER :: cosmo_lat, cosmo_lon1166 REAL(dp), DIMENSION(:,:,:), POINTER :: cosmo_h1167 1168 1178 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode) 1179 1180 TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid 1181 TYPE(grid_definition), INTENT(INOUT), TARGET :: grid, intermediate_grid 1182 CHARACTER, INTENT(IN) :: kind 1183 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 1184 1185 REAL(wp), DIMENSION(:), POINTER :: cosmo_lat, cosmo_lon 1186 REAL(wp), DIMENSION(:,:,:), POINTER :: cosmo_h 1187 1188 LOGICAL :: setup_volumetric 1169 1189 1170 1190 !------------------------------------------------------------------------------ … … 1172 1192 !------------------------------------------------------------------------------ 1173 1193 ! 1174 !-- 1175 1176 1177 ! 1178 !-- 1179 1180 1181 cosmo_lat => cosmo_grid %lat1182 cosmo_lon => cosmo_grid %lon1183 cosmo_h => cosmo_grid %hfl1184 ! 1185 !-- 1186 1187 1188 cosmo_lat => cosmo_grid %lat1189 cosmo_lon => cosmo_grid %lon1190 cosmo_h => cosmo_grid %hhl1191 ! 1192 !-- 1193 1194 1195 cosmo_lat => cosmo_grid %lat1196 cosmo_lon => cosmo_grid %lonu1197 cosmo_h => cosmo_grid %hfl1198 1199 ! 1200 !-- 1201 1202 1203 cosmo_lat => cosmo_grid %latv1204 cosmo_lon => cosmo_grid %lon1205 cosmo_h => cosmo_grid %hfl1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 intermediate_grid % clat, intermediate_grid %clon, &1216 intermediate_grid % ii, intermediate_grid %jj)1217 1218 1219 intermediate_grid % clat, intermediate_grid %clon, &1220 intermediate_grid % ii, intermediate_grid %jj, &1221 intermediate_grid %w_horiz)1194 !-- Select horizontal coordinates according to kind of points (s/w, u, v) 1195 SELECT CASE(kind) 1196 1197 ! 1198 !-- scalars 1199 CASE('s') 1200 1201 cosmo_lat => cosmo_grid%lat 1202 cosmo_lon => cosmo_grid%lon 1203 cosmo_h => cosmo_grid%hfl 1204 ! 1205 !-- vertical velocity 1206 CASE('w') 1207 1208 cosmo_lat => cosmo_grid%lat 1209 cosmo_lon => cosmo_grid%lon 1210 cosmo_h => cosmo_grid%hhl 1211 ! 1212 !-- x velocity 1213 CASE('u') 1214 1215 cosmo_lat => cosmo_grid%lat 1216 cosmo_lon => cosmo_grid%lonu 1217 cosmo_h => cosmo_grid%hfl 1218 1219 ! 1220 !-- y velocity 1221 CASE('v') 1222 1223 cosmo_lat => cosmo_grid%latv 1224 cosmo_lon => cosmo_grid%lon 1225 cosmo_h => cosmo_grid%hfl 1226 1227 CASE DEFAULT 1228 1229 message = "Interpolation quantity '" // kind // "' is not supported." 1230 CALL inifor_abort('setup_interpolation', message) 1231 1232 END SELECT 1233 1234 CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon, & 1235 intermediate_grid%clat, intermediate_grid%clon, & 1236 intermediate_grid%ii, intermediate_grid%jj) 1237 1238 CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, & 1239 intermediate_grid%clat, intermediate_grid%clon, & 1240 intermediate_grid%ii, intermediate_grid%jj, & 1241 intermediate_grid%w_horiz) 1222 1242 1223 1243 !------------------------------------------------------------------------------ … … 1226 1246 1227 1247 ! 1228 !-- 1229 !-- 1230 !-- 1231 !-- 1232 !-- 1233 1234 1235 1236 1237 1238 1239 ALLOCATE( intermediate_grid % h(0:intermediate_grid %nx, &1240 0:intermediate_grid %ny, &1241 0:intermediate_grid %nz) )1242 intermediate_grid %h(:,:,:) = - EARTH_RADIUS1243 1244 ! 1245 !-- 1246 !-- 1247 CALL interpolate_2d(cosmo_h, intermediate_grid %h, intermediate_grid)1248 1249 1250 1251 1248 !-- If profile initialization is chosen, we--somewhat counterintuitively-- 1249 !-- don't need to compute vertical interpolation weights. At least, we 1250 !-- don't need them on the intermediate grid, which fills the entire PALM 1251 !-- domain volume. Instead we need vertical weights for the intermediate 1252 !-- profile grids, which get computed in setup_averaging(). 1253 setup_volumetric = .TRUE. 1254 IF (PRESENT(ic_mode)) THEN 1255 IF (TRIM(ic_mode) == 'profile') setup_volumetric = .FALSE. 1256 ENDIF 1257 1258 IF (setup_volumetric) THEN 1259 ALLOCATE( intermediate_grid%h(0:intermediate_grid%nx, & 1260 0:intermediate_grid%ny, & 1261 0:intermediate_grid%nz) ) 1262 intermediate_grid%h(:,:,:) = - EARTH_RADIUS 1263 1264 ! 1265 !-- For w points, use hhl, for scalars use hfl 1266 !-- compute the full heights for the intermediate grids 1267 CALL interpolate_2d(cosmo_h, intermediate_grid%h, intermediate_grid) 1268 CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid) 1269 ENDIF 1270 1271 END SUBROUTINE setup_interpolation 1252 1272 1253 1273 !------------------------------------------------------------------------------! … … 1281 1301 !> grid : Grid variable to be initialized. 1282 1302 !------------------------------------------------------------------------------! 1283 1284 1285 1286 1287 1288 REAL(dp), INTENT(IN) :: xmin, xmax, ymin, ymax1289 REAL(dp), INTENT(IN) :: x0, y0, z01290 REAL(dp), INTENT(IN), TARGET, OPTIONAL :: z(:)1291 REAL(dp), INTENT(IN), TARGET, OPTIONAL :: zw(:)1292 1293 1294 grid %nx = nx1295 grid %ny = ny1296 grid %nz = nz1297 1298 grid %lx = xmax - xmin1299 grid %ly = ymax - ymin1300 1301 grid %x0 = x01302 grid %y0 = y01303 grid %z0 = z01304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 ALLOCATE( grid %x(0:nx) )1315 CALL linspace(xmin, xmax, grid %x)1316 1317 ALLOCATE( grid %y(0:ny) )1318 CALL linspace(ymin, ymax, grid %y)1319 1320 grid %z => z1321 1322 ! 1323 !-- 1324 1325 ALLOCATE( grid %kk(0:nx, 0:ny, 1:nz, 2) )1326 grid %kk(:,:,:,:) = -11327 1328 ALLOCATE( grid %w_verti(0:nx, 0:ny, 1:nz, 2) )1329 grid % w_verti(:,:,:,:) = 0.0_dp1330 1331 1332 1333 1334 ALLOCATE( grid %x(0:nx) )1335 CALL linspace(xmin, xmax, grid %x)1336 1337 ALLOCATE( grid %y(0:ny) )1338 CALL linspace(ymin, ymax, grid %y)1339 1340 ALLOCATE( grid % clon(0:nx, 0:ny), grid %clat(0:nx, 0:ny) )1341 1342 1343 phir = project( grid %y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude1344 lamr = project( grid %x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude1345 1346 phi = grid %clat, &1347 lam = grid %clon, &1348 1349 1350 1351 ! 1352 !-- 1353 ALLOCATE( grid %ii(0:nx, 0:ny, 4), &1354 grid %jj(0:nx, 0:ny, 4) )1355 grid %ii(:,:,:) = -11356 grid %jj(:,:,:) = -11357 1358 ALLOCATE( grid %w_horiz(0:nx, 0:ny, 4) )1359 grid % w_horiz(:,:,:) = 0.0_dp1360 1361 ! 1362 !-- 1363 !-- 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 grid %name(1) = 'x and lon'1377 grid %name(2) = 'y and lat'1378 grid %name(3) = 'z'1379 1380 ! 1381 !-- 1382 !-- 1383 ALLOCATE( grid % x(0:nx), grid %y(0:ny) )1384 ALLOCATE( grid % xu(1:nx), grid %yv(1:ny) )1385 CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid %x)1386 CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid %y)1387 grid %z => z1388 CALL linspace(xmin + dx, xmax - dx, grid %xu)1389 CALL linspace(ymin + dy, ymax - dy, grid %yv)1390 grid %zw => zw1391 1392 grid %depths => depths1393 1394 ! 1395 !-- 1396 1397 ALLOCATE( grid %kk(0:nx, 0:ny, 1:nz, 2) )1398 grid %kk(:,:,:,:) = -11399 1400 ALLOCATE( grid %w_verti(0:nx, 0:ny, 1:nz, 2) )1401 grid % w_verti(:,:,:,:) = 0.0_dp1402 1403 1404 1405 1406 grid %name(1) = 'x and lon'1407 grid %name(2) = 'y and lat'1408 grid %name(3) = 'interpolated hhl or hfl'1409 1410 ! 1411 !-- 1412 !-- 1413 ALLOCATE( grid % x(0:nx), grid %y(0:ny) )1414 ALLOCATE( grid % xu(1:nx), grid %yv(1:ny) )1415 CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid %x)1416 CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid %y)1417 CALL linspace(xmin + dx, xmax - dx, grid %xu)1418 CALL linspace(ymin + dy, ymax - dy, grid %yv)1419 1420 grid %depths => depths1421 1422 ! 1423 !-- 1424 ALLOCATE( grid % clon(0:nx, 0:ny), grid %clat(0:nx, 0:ny) )1425 ALLOCATE( grid % clonu(1:nx, 0:ny), grid %clatu(1:nx, 0:ny) )1426 ALLOCATE( grid % clonv(0:nx, 1:ny), grid %clatv(0:nx, 1:ny) )1427 1428 ! 1429 !-- 1430 !-- 1431 1432 phir = project( grid %y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude1433 lamr = project( grid %x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude1434 1435 phi = grid %clat, &1436 lam = grid %clon, &1437 1438 1439 1440 ! 1441 !-- 1442 1443 phir = project( grid %y, y0, EARTH_RADIUS ), & ! = plate-carree latitude1444 lamr = project( grid %xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude1445 1446 phi = grid %clatu, &1447 lam = grid %clonu, &1448 1449 1450 1451 ! 1452 !-- 1453 1454 phir = project( grid %yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude1455 lamr = project( grid %x, x0, EARTH_RADIUS ), & ! = plate-carree longitude1456 1457 phi = grid %clatv, &1458 lam = grid %clonv, &1459 1460 1461 1462 ! 1463 !-- 1464 ALLOCATE( grid %ii(0:nx, 0:ny, 4), &1465 grid %jj(0:nx, 0:ny, 4) )1466 grid %ii(:,:,:) = -11467 grid %jj(:,:,:) = -11468 1469 ALLOCATE( grid %w_horiz(0:nx, 0:ny, 4) )1470 grid % w_horiz(:,:,:) = 0.0_dp1471 1472 1473 grid %name(1) = 'rlon' ! of COMSO-DE cell centres (scalars)1474 grid %name(2) = 'rlat' ! of COMSO-DE cell centres (scalars)1475 grid %name(3) = 'height'1476 1477 ALLOCATE( grid % lon(0:nx), grid %lat(0:ny) )1478 ALLOCATE( grid % lonu(0:nx), grid %latv(0:ny) )1479 1480 CALL linspace(xmin, xmax, grid %lon)1481 CALL linspace(ymin, ymax, grid %lat)1482 grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid %nx)1483 grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid %ny)1484 1485 ! 1486 !-- 1487 !-- 1488 grid %hhl => hhl1489 grid %hfl => hfl1490 grid %depths => depths1491 1492 1493 1494 1495 1496 1497 1498 1303 SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, & 1304 x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode) 1305 CHARACTER(LEN=*), INTENT(IN) :: kind 1306 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode 1307 INTEGER, INTENT(IN) :: nx, ny, nz 1308 REAL(wp), INTENT(IN) :: xmin, xmax, ymin, ymax 1309 REAL(wp), INTENT(IN) :: x0, y0, z0 1310 REAL(wp), INTENT(IN), TARGET, OPTIONAL :: z(:) 1311 REAL(wp), INTENT(IN), TARGET, OPTIONAL :: zw(:) 1312 TYPE(grid_definition), INTENT(INOUT) :: grid 1313 1314 grid%nx = nx 1315 grid%ny = ny 1316 grid%nz = nz 1317 1318 grid%lx = xmax - xmin 1319 grid%ly = ymax - ymin 1320 1321 grid%x0 = x0 1322 grid%y0 = y0 1323 grid%z0 = z0 1324 1325 SELECT CASE( TRIM(kind) ) 1326 1327 CASE('boundary') 1328 1329 IF (.NOT.PRESENT(z)) THEN 1330 message = "z has not been passed but is required for 'boundary' grids" 1331 CALL inifor_abort('init_grid_definition', message) 1332 ENDIF 1333 1334 ALLOCATE( grid%x(0:nx) ) 1335 CALL linspace(xmin, xmax, grid%x) 1336 1337 ALLOCATE( grid%y(0:ny) ) 1338 CALL linspace(ymin, ymax, grid%y) 1339 1340 grid%z => z 1341 1342 ! 1343 !-- Allocate neighbour indices and weights 1344 IF (TRIM(ic_mode) .NE. 'profile') THEN 1345 ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) ) 1346 grid%kk(:,:,:,:) = -1 1347 1348 ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) ) 1349 grid%w_verti(:,:,:,:) = 0.0_wp 1350 ENDIF 1351 1352 CASE('boundary intermediate') 1353 1354 ALLOCATE( grid%x(0:nx) ) 1355 CALL linspace(xmin, xmax, grid%x) 1356 1357 ALLOCATE( grid%y(0:ny) ) 1358 CALL linspace(ymin, ymax, grid%y) 1359 1360 ALLOCATE( grid%clon(0:nx, 0:ny), grid%clat(0:nx, 0:ny) ) 1361 1362 CALL rotate_to_cosmo( & 1363 phir = project( grid%y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude 1364 lamr = project( grid%x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude 1365 phip = phi_cn, lamp = lambda_cn, & 1366 phi = grid%clat, & 1367 lam = grid%clon, & 1368 gam = gam & 1369 ) 1370 1371 ! 1372 !-- Allocate neighbour indices and weights 1373 ALLOCATE( grid%ii(0:nx, 0:ny, 4), & 1374 grid%jj(0:nx, 0:ny, 4) ) 1375 grid%ii(:,:,:) = -1 1376 grid%jj(:,:,:) = -1 1377 1378 ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) ) 1379 grid%w_horiz(:,:,:) = 0.0_wp 1380 1381 ! 1382 !-- This mode initializes a Cartesian PALM-4U grid and adds the 1383 !-- corresponding latitudes and longitudes of the rotated pole grid. 1384 CASE('palm') 1385 1386 IF (.NOT.PRESENT(z)) THEN 1387 message = "z has not been passed but is required for 'palm' grids" 1388 CALL inifor_abort('init_grid_definition', message) 1389 ENDIF 1390 1391 IF (.NOT.PRESENT(zw)) THEN 1392 message = "zw has not been passed but is required for 'palm' grids" 1393 CALL inifor_abort('init_grid_definition', message) 1394 ENDIF 1395 1396 grid%name(1) = 'x and lon' 1397 grid%name(2) = 'y and lat' 1398 grid%name(3) = 'z' 1399 1400 ! 1401 !-- TODO: Remove use of global dx, dy, dz variables. Consider 1402 !-- TODO: associating global x,y, and z arrays. 1403 ALLOCATE( grid%x(0:nx), grid%y(0:ny) ) 1404 ALLOCATE( grid%xu(1:nx), grid%yv(1:ny) ) 1405 CALL linspace(xmin + 0.5_wp* dx, xmax - 0.5_wp* dx, grid%x) 1406 CALL linspace(ymin + 0.5_wp* dy, ymax - 0.5_wp* dy, grid%y) 1407 grid%z => z 1408 CALL linspace(xmin + dx, xmax - dx, grid%xu) 1409 CALL linspace(ymin + dy, ymax - dy, grid%yv) 1410 grid%zw => zw 1411 1412 grid%depths => depths 1413 1414 ! 1415 !-- Allocate neighbour indices and weights 1416 IF (TRIM(ic_mode) .NE. 'profile') THEN 1417 ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) ) 1418 grid%kk(:,:,:,:) = -1 1419 1420 ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) ) 1421 grid%w_verti(:,:,:,:) = 0.0_wp 1422 ENDIF 1423 1424 CASE('palm intermediate') 1425 1426 grid%name(1) = 'x and lon' 1427 grid%name(2) = 'y and lat' 1428 grid%name(3) = 'interpolated hhl or hfl' 1429 1430 ! 1431 !-- TODO: Remove use of global dx, dy, dz variables. Consider 1432 !-- TODO: associating global x,y, and z arrays. 1433 ALLOCATE( grid%x(0:nx), grid%y(0:ny) ) 1434 ALLOCATE( grid%xu(1:nx), grid%yv(1:ny) ) 1435 CALL linspace(xmin + 0.5_wp*dx, xmax - 0.5_wp*dx, grid%x) 1436 CALL linspace(ymin + 0.5_wp*dy, ymax - 0.5_wp*dy, grid%y) 1437 CALL linspace(xmin + dx, xmax - dx, grid%xu) 1438 CALL linspace(ymin + dy, ymax - dy, grid%yv) 1439 1440 grid%depths => depths 1441 1442 ! 1443 !-- Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude 1444 ALLOCATE( grid%clon(0:nx, 0:ny), grid%clat(0:nx, 0:ny) ) 1445 ALLOCATE( grid%clonu(1:nx, 0:ny), grid%clatu(1:nx, 0:ny) ) 1446 ALLOCATE( grid%clonv(0:nx, 1:ny), grid%clatv(0:nx, 1:ny) ) 1447 1448 ! 1449 !-- Compute rotated-pole coordinates of... 1450 !-- ... PALM-4U centres 1451 CALL rotate_to_cosmo( & 1452 phir = project( grid%y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude 1453 lamr = project( grid%x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude 1454 phip = phi_cn, lamp = lambda_cn, & 1455 phi = grid%clat, & 1456 lam = grid%clon, & 1457 gam = gam & 1458 ) 1459 1460 ! 1461 !-- ... PALM-4U u winds 1462 CALL rotate_to_cosmo( & 1463 phir = project( grid%y, y0, EARTH_RADIUS ), & ! = plate-carree latitude 1464 lamr = project( grid%xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude 1465 phip = phi_cn, lamp = lambda_cn, & 1466 phi = grid%clatu, & 1467 lam = grid%clonu, & 1468 gam = gam & 1469 ) 1470 1471 ! 1472 !-- ... PALM-4U v winds 1473 CALL rotate_to_cosmo( & 1474 phir = project( grid%yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude 1475 lamr = project( grid%x, x0, EARTH_RADIUS ), & ! = plate-carree longitude 1476 phip = phi_cn, lamp = lambda_cn, & 1477 phi = grid%clatv, & 1478 lam = grid%clonv, & 1479 gam = gam & 1480 ) 1481 1482 ! 1483 !-- Allocate neighbour indices and weights 1484 ALLOCATE( grid%ii(0:nx, 0:ny, 4), & 1485 grid%jj(0:nx, 0:ny, 4) ) 1486 grid%ii(:,:,:) = -1 1487 grid%jj(:,:,:) = -1 1488 1489 ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) ) 1490 grid%w_horiz(:,:,:) = 0.0_wp 1491 1492 CASE('cosmo-de') 1493 grid%name(1) = 'rlon' ! of COMSO-DE cell centres (scalars) 1494 grid%name(2) = 'rlat' ! of COMSO-DE cell centres (scalars) 1495 grid%name(3) = 'height' 1496 1497 ALLOCATE( grid%lon(0:nx), grid%lat(0:ny) ) 1498 ALLOCATE( grid%lonu(0:nx), grid%latv(0:ny) ) 1499 1500 CALL linspace(xmin, xmax, grid%lon) 1501 CALL linspace(ymin, ymax, grid%lat) 1502 grid%lonu(:) = grid%lon + 0.5_wp * (grid%lx / grid%nx) 1503 grid%latv(:) = grid%lat + 0.5_wp * (grid%ly / grid%ny) 1504 1505 ! 1506 !-- Point to heights of half levels (hhl) and compute heights of full 1507 !-- levels (hfl) as arithmetic averages 1508 grid%hhl => hhl 1509 grid%hfl => hfl 1510 grid%depths => depths 1511 1512 CASE DEFAULT 1513 message = "Grid kind '" // TRIM(kind) // "' is not recognized." 1514 CALL inifor_abort('init_grid_definition', message) 1515 1516 END SELECT 1517 1518 END SUBROUTINE init_grid_definition 1499 1519 1500 1520 … … 1527 1547 !> avg_grid : averagin grid to be initialized 1528 1548 !------------------------------------------------------------------------------! 1529 1530 1531 1532 1533 1534 REAL(dp), INTENT(IN) :: x, y, z01535 REAL(dp), INTENT(IN), TARGET :: z(:)1536 REAL(dp), INTENT(IN) :: lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad]1537 REAL(dp), INTENT(IN) :: lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad]1538 REAL(dp), INTENT(IN) :: latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]1539 REAL(dp), INTENT(IN) :: latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]1540 1541 1542 1543 1544 1545 1546 ALLOCATE( avg_grid %x(1) )1547 ALLOCATE( avg_grid %y(1) )1548 avg_grid %x(1) = x1549 avg_grid %y(1) = y1550 avg_grid %z => z1551 avg_grid %z0 = z01552 1553 avg_grid %nz = SIZE(z, 1)1554 1555 ALLOCATE( avg_grid %lon(2) )1556 ALLOCATE( avg_grid %lat(2) )1557 avg_grid %lon(1:2) = (/lonmin, lonmax/)1558 avg_grid %lat(1:2) = (/latmin, latmax/)1559 1560 avg_grid %kind = TRIM(kind)1561 avg_grid %name(1) = TRIM(name)1562 1563 ! 1564 !-- 1565 !-- given by avg_grid % clon, %clat1566 1567 1568 ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid %nz, 2) )1569 ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid %nz, 2) )1570 ! 1571 !-- 1572 SELECT CASE(avg_grid %kind)1549 SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0, & 1550 lonmin, lonmax, latmin, latmax, kind, name) 1551 1552 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1553 TYPE(grid_definition), INTENT(IN) :: cosmo_grid 1554 REAL(wp), INTENT(IN) :: x, y, z0 1555 REAL(wp), INTENT(IN), TARGET :: z(:) 1556 REAL(wp), INTENT(IN) :: lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad] 1557 REAL(wp), INTENT(IN) :: lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad] 1558 REAL(wp), INTENT(IN) :: latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad] 1559 REAL(wp), INTENT(IN) :: latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad] 1560 1561 CHARACTER(LEN=*), INTENT(IN) :: kind 1562 CHARACTER(LEN=*), INTENT(IN) :: name 1563 1564 LOGICAL :: level_based_averaging 1565 1566 ALLOCATE( avg_grid%x(1) ) 1567 ALLOCATE( avg_grid%y(1) ) 1568 avg_grid%x(1) = x 1569 avg_grid%y(1) = y 1570 avg_grid%z => z 1571 avg_grid%z0 = z0 1572 1573 avg_grid%nz = SIZE(z, 1) 1574 1575 ALLOCATE( avg_grid%lon(2) ) 1576 ALLOCATE( avg_grid%lat(2) ) 1577 avg_grid%lon(1:2) = (/lonmin, lonmax/) 1578 avg_grid%lat(1:2) = (/latmin, latmax/) 1579 1580 avg_grid%kind = TRIM(kind) 1581 avg_grid%name(1) = TRIM(name) 1582 1583 ! 1584 !-- Find and store COSMO columns that fall into the coordinate range 1585 !-- given by avg_grid%clon, %clat 1586 CALL get_cosmo_averaging_region(avg_grid, cosmo_grid) 1587 1588 ALLOCATE (avg_grid%kkk(avg_grid%n_columns, avg_grid%nz, 2) ) 1589 ALLOCATE (avg_grid%w(avg_grid%n_columns, avg_grid%nz, 2) ) 1590 ! 1591 !-- Compute average COSMO levels in the averaging region 1592 SELECT CASE(avg_grid%kind) 1573 1593 1574 1594 CASE('scalar', 'u', 'v') 1575 avg_grid % cosmo_h => cosmo_grid %hfl1595 avg_grid%cosmo_h => cosmo_grid%hfl 1576 1596 1577 1597 CASE('w') 1578 avg_grid % cosmo_h => cosmo_grid %hhl1598 avg_grid%cosmo_h => cosmo_grid%hhl 1579 1599 1580 1600 CASE DEFAULT 1581 message = "Averaging grid kind '" // TRIM(avg_grid %kind) // &1601 message = "Averaging grid kind '" // TRIM(avg_grid%kind) // & 1582 1602 "' is not supported. Use 'scalar', 'u', or 'v'." 1583 1603 CALL inifor_abort('get_cosmo_averaging_region', message) 1584 1604 1585 END SELECT 1586 1587 ! 1588 !-- For level-besed averaging, compute average heights 1589 !level_based_averaging = ( TRIM(mode) == 'level' ) 1590 level_based_averaging = ( TRIM(cfg % averaging_mode) == 'level' ) 1591 IF (level_based_averaging) THEN 1592 ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) ) 1605 END SELECT 1606 1607 ! 1608 !-- For level-besed averaging, compute average heights 1609 level_based_averaging = ( TRIM(cfg%averaging_mode) == 'level' ) 1610 IF (level_based_averaging) THEN 1611 ALLOCATE(avg_grid%h(1,1,SIZE(avg_grid%cosmo_h, 3)) ) 1593 1612 1594 CALL average_2d(avg_grid % cosmo_h, avg_grid %h(1,1,:), &1595 avg_grid % iii, avg_grid %jjj)1596 1597 1598 1599 ! 1600 !-- 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 REAL(dp), DIMENSION(:), POINTER :: cosmo_lon, cosmo_lat1613 REAL(dp) :: dlon, dlat1614 1615 1616 1617 1618 SELECT CASE( TRIM(avg_grid %kind) )1613 CALL average_2d(avg_grid%cosmo_h, avg_grid%h(1,1,:), & 1614 avg_grid%iii, avg_grid%jjj) 1615 1616 ENDIF 1617 1618 ! 1619 !-- Compute vertical weights and neighbours 1620 CALL find_vertical_neighbours_and_weights_average( & 1621 avg_grid, level_based_averaging & 1622 ) 1623 1624 END SUBROUTINE init_averaging_grid 1625 1626 1627 SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid) 1628 TYPE(grid_definition), INTENT(INOUT) :: avg_grid 1629 TYPE(grid_definition), TARGET, INTENT(IN) :: cosmo_grid 1630 1631 REAL(wp), DIMENSION(:), POINTER :: cosmo_lon, cosmo_lat 1632 REAL(wp) :: dlon, dlat 1633 1634 INTEGER :: i, j, imin, imax, jmin, jmax, l, nx, ny 1635 1636 1637 SELECT CASE( TRIM(avg_grid%kind) ) 1619 1638 1620 1639 CASE('scalar', 'w') 1621 cosmo_lon => cosmo_grid %lon1622 cosmo_lat => cosmo_grid %lat1640 cosmo_lon => cosmo_grid%lon 1641 cosmo_lat => cosmo_grid%lat 1623 1642 1624 1643 CASE('u') 1625 cosmo_lon => cosmo_grid %lonu1626 cosmo_lat => cosmo_grid %lat1644 cosmo_lon => cosmo_grid%lonu 1645 cosmo_lat => cosmo_grid%lat 1627 1646 1628 1647 CASE('v') 1629 cosmo_lon => cosmo_grid %lon1630 cosmo_lat => cosmo_grid %latv1648 cosmo_lon => cosmo_grid%lon 1649 cosmo_lat => cosmo_grid%latv 1631 1650 1632 1651 CASE DEFAULT 1633 message = "Averaging grid kind '" // TRIM(avg_grid %kind) // &1652 message = "Averaging grid kind '" // TRIM(avg_grid%kind) // & 1634 1653 "' is not supported. Use 'scalar', 'u', or 'v'." 1635 1654 CALL inifor_abort('get_cosmo_averaging_region', message) 1636 1655 1637 1638 1639 ! 1640 !-- 1641 1642 1643 1644 imin = FLOOR ( (avg_grid %lon(1) - cosmo_lon(0)) / dlon )1645 imax = CEILING( (avg_grid %lon(2) - cosmo_lon(0)) / dlon )1646 1647 jmin = FLOOR ( (avg_grid %lat(1) - cosmo_lat(0)) / dlat )1648 jmax = CEILING( (avg_grid %lat(2) - cosmo_lat(0)) / dlat )1649 1650 message = "Grid " // TRIM(avg_grid %name(1)) // " averages over " // &1651 1652 1653 1654 1655 1656 1657 1658 avg_grid %n_columns = nx * ny1659 1660 ALLOCATE( avg_grid % iii(avg_grid %n_columns), &1661 avg_grid % jjj(avg_grid %n_columns) )1662 1663 1664 DOj = jmin, jmax1665 DOi = imin, imax1666 1667 avg_grid %iii(l) = i1668 avg_grid %jjj(l) = j1669 1670 1671 1672 1656 END SELECT 1657 1658 ! 1659 !-- FIXME: make dlon, dlat parameters of the grid_defintion type 1660 dlon = cosmo_lon(1) - cosmo_lon(0) 1661 dlat = cosmo_lat(1) - cosmo_lat(0) 1662 1663 imin = FLOOR ( (avg_grid%lon(1) - cosmo_lon(0)) / dlon ) 1664 imax = CEILING( (avg_grid%lon(2) - cosmo_lon(0)) / dlon ) 1665 1666 jmin = FLOOR ( (avg_grid%lat(1) - cosmo_lat(0)) / dlat ) 1667 jmax = CEILING( (avg_grid%lat(2) - cosmo_lat(0)) / dlat ) 1668 1669 message = "Grid " // TRIM(avg_grid%name(1)) // " averages over " // & 1670 TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) // & 1671 " and " // & 1672 TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax)) 1673 CALL report( 'get_cosmo_averaging_region', message ) 1674 1675 nx = imax - imin + 1 1676 ny = jmax - jmin + 1 1677 avg_grid%n_columns = nx * ny 1678 1679 ALLOCATE( avg_grid%iii(avg_grid%n_columns), & 1680 avg_grid%jjj(avg_grid%n_columns) ) 1681 1682 l = 0 1683 DO j = jmin, jmax 1684 DO i = imin, imax 1685 l = l + 1 1686 avg_grid%iii(l) = i 1687 avg_grid%jjj(l) = j 1688 ENDDO 1689 ENDDO 1690 1691 END SUBROUTINE get_cosmo_averaging_region 1673 1692 1674 1693 … … 1683 1702 !> 'modpoints'. 1684 1703 !------------------------------------------------------------------------------! 1685 1686 1687 1688 1689 REAL(dp), DIMENSION(:), INTENT(INOUT) :: z, dz, dz_stretch_factor_array1690 REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_level_start, dz_stretch_level_end1691 REAL(dp), INTENT(IN) :: dz_max, dz_stretch_factor, dz_stretch_level1692 1693 1694 1695 1696 REAL(dp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end1697 REAL(dp) :: dz_level_end, dz_stretched1698 1699 1700 1701 1702 1704 SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, & 1705 dz_stretch_level_start, dz_stretch_level_end, & 1706 dz_stretch_factor_array) 1707 1708 REAL(wp), DIMENSION(:), INTENT(INOUT) :: z, dz, dz_stretch_factor_array 1709 REAL(wp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_level_start, dz_stretch_level_end 1710 REAL(wp), INTENT(IN) :: dz_max, dz_stretch_factor, dz_stretch_level 1711 1712 INTEGER :: number_stretch_level_start !< number of user-specified start levels for stretching 1713 INTEGER :: number_stretch_level_end !< number of user-specified end levels for stretching 1714 1715 REAL(wp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end 1716 REAL(wp) :: dz_level_end, dz_stretched 1717 1718 INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 1719 INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 1720 INTEGER :: dz_stretch_level_index = 0 1721 INTEGER :: k, n, number_dz 1703 1722 1704 1723 ! 1705 1724 !-- Compute height of u-levels from constant grid length and dz stretch factors 1706 IF ( dz(1) == -1.0_dp ) THEN1707 1708 1709 ELSEIF ( dz(1) <= 0.0_dp ) THEN1710 1711 1712 1725 IF ( dz(1) == -1.0_wp ) THEN 1726 message = 'missing dz' 1727 CALL inifor_abort( 'stretched_z', message) 1728 ELSEIF ( dz(1) <= 0.0_wp ) THEN 1729 WRITE( message, * ) 'dz=', dz(1),' <= 0.0' 1730 CALL inifor_abort( 'stretched_z', message) 1731 ENDIF 1713 1732 1714 1733 ! 1715 1734 !-- Initialize dz_stretch_level_start with the value of dz_stretch_level 1716 1735 !-- if it was set by the user 1717 IF ( dz_stretch_level /= -9999999.9_dp )THEN1718 1719 1736 IF ( dz_stretch_level /= -9999999.9_wp ) THEN 1737 dz_stretch_level_start(1) = dz_stretch_level 1738 ENDIF 1720 1739 1721 1740 ! … … 1727 1746 !-- is used (Attention: The user is not allowed to specify a dz value equal 1728 1747 !-- to the default of dz_max = 999.0). 1729 number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max )1730 1731 -9999999.9_dp )1732 1733 9999999.9_dp )1748 number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max ) 1749 number_stretch_level_start = COUNT( dz_stretch_level_start /= & 1750 -9999999.9_wp ) 1751 number_stretch_level_end = COUNT( dz_stretch_level_end /= & 1752 9999999.9_wp ) 1734 1753 1735 1754 ! 1736 1755 !-- The number of specified end levels +1 has to be the same than the number 1737 1756 !-- of specified dz values 1738 IF ( number_dz /= number_stretch_level_end + 1 )THEN1739 1740 1741 1742 1743 1744 1745 1757 IF ( number_dz /= number_stretch_level_end + 1 ) THEN 1758 WRITE( message, * ) 'The number of values for dz = ', & 1759 number_dz, 'has to be the same than ', & 1760 'the number of values for ', & 1761 'dz_stretch_level_end + 1 = ', & 1762 number_stretch_level_end+1 1763 CALL inifor_abort( 'stretched_z', message) 1764 ENDIF 1746 1765 1747 1766 ! 1748 !-- 1749 !-- 1750 1751 number_dz /= number_stretch_level_start )THEN1752 1753 1754 1755 1756 1757 1758 1767 !-- The number of specified start levels has to be the same or one less than 1768 !-- the number of specified dz values 1769 IF ( number_dz /= number_stretch_level_start + 1 .AND. & 1770 number_dz /= number_stretch_level_start ) THEN 1771 WRITE( message, * ) 'The number of values for dz = ', & 1772 number_dz, 'has to be the same or one ', & 1773 'more than& the number of values for ', & 1774 'dz_stretch_level_start = ', & 1775 number_stretch_level_start 1776 CALL inifor_abort( 'stretched_z', message) 1777 ENDIF 1759 1778 1760 !-- 1761 !-- 1762 1763 number_stretch_level_start /= number_stretch_level_end )THEN1764 1765 1766 1767 1768 1769 1770 1771 1779 !-- The number of specified start levels has to be the same or one more than 1780 !-- the number of specified end levels 1781 IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. & 1782 number_stretch_level_start /= number_stretch_level_end ) THEN 1783 WRITE( message, * ) 'The number of values for ', & 1784 'dz_stretch_level_start = ', & 1785 dz_stretch_level_start, 'has to be the ',& 1786 'same or one more than& the number of ', & 1787 'values for dz_stretch_level_end = ', & 1788 number_stretch_level_end 1789 CALL inifor_abort( 'stretched_z', message) 1790 ENDIF 1772 1791 1773 1792 ! 1774 1793 !-- Initialize dz for the free atmosphere with the value of dz_max 1775 IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND. &1776 number_stretch_level_start /= 0 )THEN1777 1778 1794 IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND. & 1795 number_stretch_level_start /= 0 ) THEN 1796 dz(number_stretch_level_start+1) = dz_max 1797 ENDIF 1779 1798 1780 1799 ! … … 1782 1801 !-- atmosphere is desired (dz_stretch_level_end was not specified for the 1783 1802 !-- free atmosphere) 1784 IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 1785 dz_stretch_factor_array(number_stretch_level_start) = & 1786 dz_stretch_factor 1803 IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 1804 dz_stretch_factor_array(number_stretch_level_start) = & 1805 dz_stretch_factor 1806 ENDIF 1807 1808 !-- Allocation of arrays for stretching 1809 ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) ) 1810 1811 ! 1812 !-- The stretching region has to be large enough to allow for a smooth 1813 !-- transition between two different grid spacings 1814 DO n = 1, number_stretch_level_start 1815 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & 1816 4 * MAX( dz(n),dz(n+1) ) 1817 ENDDO 1818 1819 IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > & 1820 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN 1821 !IF ( ANY( min_dz_stretch_level_end > & 1822 ! dz_stretch_level_end ) ) THEN 1823 message = 'Each dz_stretch_level_end has to be larger ' // & 1824 'than its corresponding value for ' // & 1825 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& 1826 'to allow for smooth grid stretching' 1827 CALL inifor_abort('stretched_z', message) 1828 ENDIF 1829 1830 ! 1831 !-- Stretching must not be applied within the prandtl_layer 1832 !-- (first two grid points). For the default case dz_stretch_level_start 1833 !-- is negative. Therefore the absolut value is checked here. 1834 IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) ) THEN 1835 WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',& 1836 'larger than ', dz(1) * 1.5 1837 CALL inifor_abort( 'stretched_z', message) 1838 ENDIF 1839 1840 ! 1841 !-- The stretching has to start and end on a grid level. Therefore 1842 !-- user-specified values have to ''interpolate'' to the next lowest level 1843 IF ( number_stretch_level_start /= 0 ) THEN 1844 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & 1845 dz(1)/2.0) / dz(1) ) & 1846 * dz(1) + dz(1)/2.0 1847 ENDIF 1848 1849 IF ( number_stretch_level_start > 1 ) THEN 1850 DO n = 2, number_stretch_level_start 1851 dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & 1852 dz(n) ) * dz(n) 1853 ENDDO 1854 ENDIF 1855 1856 IF ( number_stretch_level_end /= 0 ) THEN 1857 DO n = 1, number_stretch_level_end 1858 dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & 1859 dz(n+1) ) * dz(n+1) 1860 ENDDO 1861 ENDIF 1862 1863 ! 1864 !-- Determine stretching factor if necessary 1865 IF ( number_stretch_level_end >= 1 ) THEN 1866 CALL calculate_stretching_factor( number_stretch_level_end, dz, & 1867 dz_stretch_factor_array, & 1868 dz_stretch_level_end, & 1869 dz_stretch_level_start ) 1870 ENDIF 1871 1872 z(1) = dz(1) * 0.5_wp 1873 ! 1874 dz_stretch_level_index = n 1875 dz_stretched = dz(1) 1876 DO k = 2, n 1877 1878 IF ( dz_stretch_level <= z(k-1) .AND. dz_stretched < dz_max ) THEN 1879 1880 dz_stretched = dz_stretched * dz_stretch_factor 1881 dz_stretched = MIN( dz_stretched, dz_max ) 1882 1883 IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1 1884 1787 1885 ENDIF 1788 1886 1789 !-- Allocation of arrays for stretching 1790 ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) ) 1791 1792 ! 1793 !-- The stretching region has to be large enough to allow for a smooth 1794 !-- transition between two different grid spacings 1795 DO n = 1, number_stretch_level_start 1796 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & 1797 4 * MAX( dz(n),dz(n+1) ) 1798 ENDDO 1799 1800 IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > & 1801 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN 1802 !IF ( ANY( min_dz_stretch_level_end > & 1803 ! dz_stretch_level_end ) ) THEN 1804 message = 'Each dz_stretch_level_end has to be larger ' // & 1805 'than its corresponding value for ' // & 1806 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& 1807 'to allow for smooth grid stretching' 1808 CALL inifor_abort('stretched_z', message) 1809 ENDIF 1810 1811 ! 1812 !-- Stretching must not be applied within the prandtl_layer 1813 !-- (first two grid points). For the default case dz_stretch_level_start 1814 !-- is negative. Therefore the absolut value is checked here. 1815 IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN 1816 WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',& 1817 'larger than ', dz(1) * 1.5 1818 CALL inifor_abort( 'stretched_z', message) 1819 ENDIF 1820 1821 ! 1822 !-- The stretching has to start and end on a grid level. Therefore 1823 !-- user-specified values have to ''interpolate'' to the next lowest level 1824 IF ( number_stretch_level_start /= 0 ) THEN 1825 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & 1826 dz(1)/2.0) / dz(1) ) & 1827 * dz(1) + dz(1)/2.0 1828 ENDIF 1829 1830 IF ( number_stretch_level_start > 1 ) THEN 1831 DO n = 2, number_stretch_level_start 1832 dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & 1833 dz(n) ) * dz(n) 1834 ENDDO 1835 ENDIF 1836 1837 IF ( number_stretch_level_end /= 0 ) THEN 1838 DO n = 1, number_stretch_level_end 1839 dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & 1840 dz(n+1) ) * dz(n+1) 1841 ENDDO 1842 ENDIF 1843 1844 ! 1845 !-- Determine stretching factor if necessary 1846 IF ( number_stretch_level_end >= 1 ) THEN 1847 CALL calculate_stretching_factor( number_stretch_level_end, dz, & 1848 dz_stretch_factor_array, & 1849 dz_stretch_level_end, & 1850 dz_stretch_level_start ) 1851 ENDIF 1852 1853 z(1) = dz(1) * 0.5_dp 1854 ! 1855 dz_stretch_level_index = n 1856 dz_stretched = dz(1) 1857 DO k = 2, n 1858 1859 IF ( dz_stretch_level <= z(k-1) .AND. dz_stretched < dz_max ) THEN 1860 1861 dz_stretched = dz_stretched * dz_stretch_factor 1862 dz_stretched = MIN( dz_stretched, dz_max ) 1863 1864 IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1 1865 1866 ENDIF 1867 1868 z(k) = z(k-1) + dz_stretched 1869 1870 ENDDO 1871 !-- Determine u and v height levels considering the possibility of grid 1872 !-- stretching in several heights. 1873 n = 1 1874 dz_stretch_level_start_index(:) = UBOUND(z, 1) 1875 dz_stretch_level_end_index(:) = UBOUND(z, 1) 1876 dz_stretched = dz(1) 1877 1878 !-- The default value of dz_stretch_level_start is negative, thus the first 1879 !-- condition is always true. Hence, the second condition is necessary. 1880 DO k = 2, UBOUND(z, 1) 1881 IF ( dz_stretch_level_start(n) <= z(k-1) .AND. & 1882 dz_stretch_level_start(n) /= -9999999.9_dp ) THEN 1883 dz_stretched = dz_stretched * dz_stretch_factor_array(n) 1884 1885 IF ( dz(n) > dz(n+1) ) THEN 1886 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz 1887 ELSE 1888 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz 1889 ENDIF 1890 1891 IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) ) & 1892 dz_stretch_level_start_index(n) = k-1 1893 1887 z(k) = z(k-1) + dz_stretched 1888 1889 ENDDO 1890 !-- Determine u and v height levels considering the possibility of grid 1891 !-- stretching in several heights. 1892 n = 1 1893 dz_stretch_level_start_index(:) = UBOUND(z, 1) 1894 dz_stretch_level_end_index(:) = UBOUND(z, 1) 1895 dz_stretched = dz(1) 1896 1897 !-- The default value of dz_stretch_level_start is negative, thus the first 1898 !-- condition is always true. Hence, the second condition is necessary. 1899 DO k = 2, UBOUND(z, 1) 1900 IF ( dz_stretch_level_start(n) <= z(k-1) .AND. & 1901 dz_stretch_level_start(n) /= -9999999.9_wp ) THEN 1902 dz_stretched = dz_stretched * dz_stretch_factor_array(n) 1903 1904 IF ( dz(n) > dz(n+1) ) THEN 1905 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz 1906 ELSE 1907 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz 1894 1908 ENDIF 1895 1909 1896 z(k) = z(k-1) + dz_stretched 1910 IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) ) & 1911 dz_stretch_level_start_index(n) = k-1 1897 1912 1898 ! 1899 !-- Make sure that the stretching ends exactly at dz_stretch_level_end 1900 dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) 1901 1902 IF ( dz_level_end < dz(n+1)/3.0 ) THEN 1903 z(k) = dz_stretch_level_end(n) 1904 dz_stretched = dz(n+1) 1905 dz_stretch_level_end_index(n) = k 1906 n = n + 1 1907 ENDIF 1908 ENDDO 1909 1910 DEALLOCATE( min_dz_stretch_level_end ) 1911 1912 END SUBROUTINE stretched_z 1913 ENDIF 1914 1915 z(k) = z(k-1) + dz_stretched 1916 1917 ! 1918 !-- Make sure that the stretching ends exactly at dz_stretch_level_end 1919 dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) 1920 1921 IF ( dz_level_end < dz(n+1)/3.0 ) THEN 1922 z(k) = dz_stretch_level_end(n) 1923 dz_stretched = dz(n+1) 1924 dz_stretch_level_end_index(n) = k 1925 n = n + 1 1926 ENDIF 1927 ENDDO 1928 1929 DEALLOCATE( min_dz_stretch_level_end ) 1930 1931 END SUBROUTINE stretched_z 1913 1932 1914 1933 … … 1928 1947 dz_stretch_level_start ) 1929 1948 1930 REAL( dp), DIMENSION(:), INTENT(IN) :: dz1931 REAL( dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_factor_array1932 REAL( dp), DIMENSION(:), INTENT(IN) :: dz_stretch_level_end, dz_stretch_level_start1949 REAL(wp), DIMENSION(:), INTENT(IN) :: dz 1950 REAL(wp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_factor_array 1951 REAL(wp), DIMENSION(:), INTENT(IN) :: dz_stretch_level_end, dz_stretch_level_start 1933 1952 1934 1953 INTEGER :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached … … 1938 1957 INTEGER, INTENT(IN) :: number_end !< number of user-specified end levels for stretching 1939 1958 1940 REAL( dp) :: delta_l !< absolute difference between l and l_rounded1941 REAL( dp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_21942 REAL( dp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)1943 REAL( dp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration1944 REAL( dp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)1945 REAL( dp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly1946 REAL( dp) :: numerator !< numerator of the quotient1947 REAL( dp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly1948 REAL( dp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly1959 REAL(wp) :: delta_l !< absolute difference between l and l_rounded 1960 REAL(wp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2 1961 REAL(wp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) 1962 REAL(wp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration 1963 REAL(wp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region) 1964 REAL(wp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly 1965 REAL(wp) :: numerator !< numerator of the quotient 1966 REAL(wp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly 1967 REAL(wp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly 1949 1968 1950 REAL( dp) :: dz_stretch_factor_array_2(9) = 1.08_dp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_11969 REAL(wp) :: dz_stretch_factor_array_2(9) = 1.08_wp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 1951 1970 1952 REAL( dp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors1953 REAL( dp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor1954 REAL( dp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor1971 REAL(wp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors 1972 REAL(wp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor 1973 REAL(wp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor 1955 1974 1956 1975 … … 1963 1982 delta_total_old = 1.0 1964 1983 1965 IF ( dz(n) > dz(n+1) ) THEN1966 DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )1984 IF ( dz(n) > dz(n+1) ) THEN 1985 DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 1967 1986 1968 1987 stretch_factor_1 = 1.0 - iterations * stretch_factor_interval … … 1972 1991 stretch_factor_1 - distance/dz(n) 1973 1992 1974 IF ( numerator > 0.0 ) THEN1993 IF ( numerator > 0.0 ) THEN 1975 1994 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1976 1995 l_rounded = NINT( l ) … … 1991 2010 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1992 2011 !-- equal to dz(n+1) after l_rounded grid levels. 1993 IF (delta_total_new < delta_total_old) THEN2012 IF (delta_total_new < delta_total_old) THEN 1994 2013 dz_stretch_factor_array(n) = stretch_factor_1 1995 2014 dz_stretch_factor_array_2(n) = stretch_factor_2 … … 2001 2020 ENDDO 2002 2021 2003 ELSEIF ( dz(n) < dz(n+1) ) THEN2004 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )2022 ELSEIF ( dz(n) < dz(n+1) ) THEN 2023 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) 2005 2024 2006 2025 stretch_factor_1 = 1.0 + iterations * stretch_factor_interval … … 2027 2046 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 2028 2047 !-- equal to dz(n+1) after l_rounded grid levels. 2029 IF (delta_total_new < delta_total_old) THEN2048 IF (delta_total_new < delta_total_old) THEN 2030 2049 dz_stretch_factor_array(n) = stretch_factor_1 2031 2050 dz_stretch_factor_array_2(n) = stretch_factor_2 … … 2045 2064 !-- interval. If not, print a warning for the user. 2046 2065 IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR. & 2047 dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN2048 WRITE( message, * ) 'stretch_factor_2 = ', &2066 dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN 2067 WRITE( message, * ) 'stretch_factor_2 = ', & 2049 2068 dz_stretch_factor_array_2(n), ' which is',& 2050 2069 ' responsible for exactly reaching& dz =',& … … 2067 2086 !> coordinate vector 'z' and stores it in 'zw'. 2068 2087 !------------------------------------------------------------------------------! 2069 2070 2071 REAL(dp), INTENT(IN) :: z(0:)2072 REAL(dp), INTENT(OUT) :: zw(1:)2073 2074 2075 2076 DOk = 1, UBOUND(zw, 1)2077 zw(k) = 0.5_dp * (z(k-1) + z(k))2078 2079 2080 2088 SUBROUTINE midpoints(z, zw) 2089 2090 REAL(wp), INTENT(IN) :: z(0:) 2091 REAL(wp), INTENT(OUT) :: zw(1:) 2092 2093 INTEGER :: k 2094 2095 DO k = 1, UBOUND(zw, 1) 2096 zw(k) = 0.5_wp * (z(k-1) + z(k)) 2097 ENDDO 2098 2099 END SUBROUTINE midpoints 2081 2100 2082 2101 !------------------------------------------------------------------------------! … … 2085 2104 !> Defines INFOR's IO groups. 2086 2105 !------------------------------------------------------------------------------! 2087 2088 2089 2090 2091 2092 2093 2094 ! 2095 !-- 2096 2097 2098 2099 2100 2101 2102 2103 ! 2104 !-- 2105 2106 2107 2108 2109 2110 2111 2112 ! 2113 !-- 2114 !-- 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127 ! 2128 !-- 2129 2130 2131 2132 2133 2134 2135 2136 ! 2137 !-- 2138 2139 2140 2141 2142 2143 2144 2145 2106 SUBROUTINE setup_io_groups() 2107 2108 INTEGER :: ngroups 2109 2110 ngroups = 16 2111 ALLOCATE( io_group_list(ngroups) ) 2112 2113 ! 2114 !-- soil temp 2115 io_group_list(1) = init_io_group( & 2116 in_files = soil_files, & 2117 out_vars = output_var_table(1:1), & 2118 in_var_list = input_var_table(1:1), & 2119 kind = 'soil-temperature' & 2120 ) 2121 2122 ! 2123 !-- soil water 2124 io_group_list(2) = init_io_group( & 2125 in_files = soil_files, & 2126 out_vars = output_var_table(2:2), & 2127 in_var_list = input_var_table(2:2), & 2128 kind = 'soil-water' & 2129 ) 2130 2131 ! 2132 !-- potential temperature, surface pressure, specific humidity including 2133 !-- nudging and subsidence, and geostrophic winds ug, vg 2134 io_group_list(3) = init_io_group( & 2135 in_files = flow_files, & 2136 out_vars = [output_var_table(56:64), & ! internal averaged density and pressure profiles 2137 output_var_table(3:8), output_var_table(9:14), & 2138 output_var_table(42:42), output_var_table(43:44), & 2139 output_var_table(49:51), output_var_table(52:54)], & 2140 in_var_list = (/input_var_table(3), input_var_table(17), & ! T, P, QV 2141 input_var_table(4) /), & 2142 kind = 'thermodynamics', & 2143 n_output_quantities = 4 & ! P, Theta, Rho, qv 2144 ) 2145 2146 ! 2147 !-- Moved to therodynamic io_group 2148 !io_group_list(4) = init_io_group( & 2149 ! in_files = flow_files, & 2150 ! out_vars = [output_var_table(9:14), output_var_table(52:54)], & 2151 ! in_var_list = input_var_table(4:4), & 2152 ! kind = 'scalar' & 2153 !) 2154 2155 ! 2156 !-- u and v velocity 2157 io_group_list(5) = init_io_group( & 2158 in_files = flow_files, & 2159 out_vars = [output_var_table(15:26), output_var_table(45:46)], & 2160 !out_vars = output_var_table(15:20), & 2161 in_var_list = input_var_table(5:6), & 2162 !in_var_list = input_var_table(5:5), & 2163 kind = 'velocities' & 2164 ) 2146 2165 2147 2166 ! 2148 !-- 2149 2150 2151 2152 2153 2154 2155 !io_group_list(6) %to_be_processed = .FALSE.2167 !-- v velocity, deprecated! 2168 !io_group_list(6) = init_io_group( & 2169 ! in_files = flow_files, & 2170 ! out_vars = output_var_table(21:26), & 2171 ! in_var_list = input_var_table(6:6), & 2172 ! kind = 'horizontal velocity' & 2173 !) 2174 !io_group_list(6)%to_be_processed = .FALSE. 2156 2175 2157 2176 ! 2158 !-- 2159 2160 2161 2162 2163 2164 2165 ! 2166 !-- 2167 2168 2169 2170 2171 2172 2173 io_group_list(8) %to_be_processed = .FALSE.2174 ! 2175 !-- 2176 2177 2178 2179 2180 2181 2182 io_group_list(9) %to_be_processed = .FALSE.2183 ! 2184 !-- 2185 2186 2187 2188 2189 2190 2191 io_group_list(10) %to_be_processed = .FALSE.2192 ! 2193 !-- 2194 2195 2196 2197 2198 2199 2200 io_group_list(11) %to_be_processed = .FALSE.2201 ! 2202 !-- 2203 2204 2205 2206 2207 2208 2209 io_group_list(12) %to_be_processed = .FALSE.2210 ! 2211 !-- 2212 2213 2214 2215 2216 2217 2218 io_group_list(13) %to_be_processed = .FALSE.2219 ! 2220 !-- 2221 2222 2223 2224 2225 2226 2227 io_group_list(14) %to_be_processed = .FALSE.2228 ! 2229 !-- 2230 2231 2232 2233 2234 2235 2236 io_group_list(15) %to_be_processed = .FALSE.2237 ! 2238 !-- 2239 2240 2241 2242 2243 2244 2245 io_group_list(16) %to_be_processed = .FALSE.2246 2247 2177 !-- w velocity and subsidence and w nudging 2178 io_group_list(7) = init_io_group( & 2179 in_files = flow_files, & 2180 out_vars = [output_var_table(27:32), output_var_table(47:48)], & 2181 in_var_list = input_var_table(7:7), & 2182 kind = 'scalar' & 2183 ) 2184 ! 2185 !-- rain 2186 io_group_list(8) = init_io_group( & 2187 in_files = soil_moisture_files, & 2188 out_vars = output_var_table(33:33), & 2189 in_var_list = input_var_table(8:8), & 2190 kind = 'accumulated' & 2191 ) 2192 io_group_list(8)%to_be_processed = .FALSE. 2193 ! 2194 !-- snow 2195 io_group_list(9) = init_io_group( & 2196 in_files = soil_moisture_files, & 2197 out_vars = output_var_table(34:34), & 2198 in_var_list = input_var_table(9:9), & 2199 kind = 'accumulated' & 2200 ) 2201 io_group_list(9)%to_be_processed = .FALSE. 2202 ! 2203 !-- graupel 2204 io_group_list(10) = init_io_group( & 2205 in_files = soil_moisture_files, & 2206 out_vars = output_var_table(35:35), & 2207 in_var_list = input_var_table(10:10), & 2208 kind = 'accumulated' & 2209 ) 2210 io_group_list(10)%to_be_processed = .FALSE. 2211 ! 2212 !-- evapotranspiration 2213 io_group_list(11) = init_io_group( & 2214 in_files = soil_moisture_files, & 2215 out_vars = output_var_table(37:37), & 2216 in_var_list = input_var_table(11:11), & 2217 kind = 'accumulated' & 2218 ) 2219 io_group_list(11)%to_be_processed = .FALSE. 2220 ! 2221 !-- 2m air temperature 2222 io_group_list(12) = init_io_group( & 2223 in_files = soil_moisture_files, & 2224 out_vars = output_var_table(36:36), & 2225 in_var_list = input_var_table(12:12), & 2226 kind = 'surface' & 2227 ) 2228 io_group_list(12)%to_be_processed = .FALSE. 2229 ! 2230 !-- incoming diffusive sw flux 2231 io_group_list(13) = init_io_group( & 2232 in_files = radiation_files, & 2233 out_vars = output_var_table(38:38), & 2234 in_var_list = input_var_table(13:13), & 2235 kind = 'running average' & 2236 ) 2237 io_group_list(13)%to_be_processed = .FALSE. 2238 ! 2239 !-- incoming direct sw flux 2240 io_group_list(14) = init_io_group( & 2241 in_files = radiation_files, & 2242 out_vars = output_var_table(39:39), & 2243 in_var_list = input_var_table(14:14), & 2244 kind = 'running average' & 2245 ) 2246 io_group_list(14)%to_be_processed = .FALSE. 2247 ! 2248 !-- sw radiation balance 2249 io_group_list(15) = init_io_group( & 2250 in_files = radiation_files, & 2251 out_vars = output_var_table(40:40), & 2252 in_var_list = input_var_table(15:15), & 2253 kind = 'running average' & 2254 ) 2255 io_group_list(15)%to_be_processed = .FALSE. 2256 ! 2257 !-- lw radiation balance 2258 io_group_list(16) = init_io_group( & 2259 in_files = radiation_files, & 2260 out_vars = output_var_table(41:41), & 2261 in_var_list = input_var_table(16:16), & 2262 kind = 'running average' & 2263 ) 2264 io_group_list(16)%to_be_processed = .FALSE. 2265 2266 END SUBROUTINE setup_io_groups 2248 2267 2249 2268 … … 2259 2278 !> on the output quantity Theta. 2260 2279 !------------------------------------------------------------------------------! 2261 2262 2263 2264 2265 2266 2267 2268 2269 2270 2271 group %nt = SIZE(in_files)2272 group %nv = SIZE(out_vars)2273 group %n_inputs = SIZE(in_var_list)2274 group %kind = TRIM(kind)2275 ! 2276 !-- 2277 !-- 2278 !-- 2279 !-- 2280 !-- 2281 !-- 2282 2283 group %n_output_quantities = n_output_quantities2284 2285 group % n_output_quantities = group %n_inputs2286 2287 2288 ALLOCATE(group % in_var_list(group %n_inputs))2289 ALLOCATE(group % in_files(group %nt))2290 ALLOCATE(group % out_vars(group %nv))2291 2292 group %in_var_list = in_var_list2293 group %in_files = in_files2294 group %out_vars = out_vars2295 group %to_be_processed = .TRUE.2296 2297 2280 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind, & 2281 n_output_quantities) RESULT(group) 2282 CHARACTER(LEN=PATH), INTENT(IN) :: in_files(:) 2283 CHARACTER(LEN=*), INTENT(IN) :: kind 2284 TYPE(nc_var), INTENT(IN) :: out_vars(:) 2285 TYPE(nc_var), INTENT(IN) :: in_var_list(:) 2286 INTEGER, OPTIONAL :: n_output_quantities 2287 2288 TYPE(io_group) :: group 2289 2290 group%nt = SIZE(in_files) 2291 group%nv = SIZE(out_vars) 2292 group%n_inputs = SIZE(in_var_list) 2293 group%kind = TRIM(kind) 2294 ! 2295 !-- For the 'thermodynamics' IO group, one quantity more than input variables 2296 !-- is needed to compute all output variables of the IO group. Concretely, in 2297 !-- preprocess() the density is computed from T,P or PP,QV in adddition to 2298 !-- the variables Theta, p, qv. In read_input_variables(), 2299 !-- n_output_quantities is used to allocate the correct number of input 2300 !-- buffers. 2301 IF ( PRESENT(n_output_quantities) ) THEN 2302 group%n_output_quantities = n_output_quantities 2303 ELSE 2304 group%n_output_quantities = group%n_inputs 2305 ENDIF 2306 2307 ALLOCATE(group%in_var_list(group%n_inputs)) 2308 ALLOCATE(group%in_files(group%nt)) 2309 ALLOCATE(group%out_vars(group%nv)) 2310 2311 group%in_var_list = in_var_list 2312 group%in_files = in_files 2313 group%out_vars = out_vars 2314 group%to_be_processed = .TRUE. 2315 2316 END FUNCTION init_io_group 2298 2317 2299 2318 … … 2303 2322 !> Deallocates all allocated variables. 2304 2323 !------------------------------------------------------------------------------! 2305 2306 2307 CALL report('fini_grids', 'Deallocating grids', cfg %debug)2308 2309 2310 2311 DEALLOCATE(palm_grid%x, palm_grid%y, palm_grid%z,&2312 palm_grid%xu, palm_grid%yv, palm_grid%zw,&2313 palm_grid%clon, palm_grid%clat,&2314 2315 2316 2317 2318 palm_intermediate%clon, palm_intermediate%clat,&2319 2320 2321 DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat,&2322 cosmo_grid%lonu, cosmo_grid%latv,&2323 2324 2325 2324 SUBROUTINE fini_grids() 2325 2326 CALL report('fini_grids', 'Deallocating grids', cfg%debug) 2327 2328 DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column) 2329 2330 DEALLOCATE(palm_grid%x, palm_grid%y, palm_grid%z, & 2331 palm_grid%xu, palm_grid%yv, palm_grid%zw, & 2332 palm_grid%clon, palm_grid%clat, & 2333 palm_grid%clonu, palm_grid%clatu) 2334 2335 DEALLOCATE(palm_intermediate%x, palm_intermediate%y, palm_intermediate%z, & 2336 palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,& 2337 palm_intermediate%clon, palm_intermediate%clat, & 2338 palm_intermediate%clonu, palm_intermediate%clatu) 2339 2340 DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat, & 2341 cosmo_grid%lonu, cosmo_grid%latv, & 2342 cosmo_grid%hfl) 2343 2344 END SUBROUTINE fini_grids 2326 2345 2327 2346 … … 2331 2350 !> Initializes the variable list. 2332 2351 !------------------------------------------------------------------------------! 2333 2334 2335 2336 2337 2338 2339 IF (TRIM(cfg %start_date) == '') THEN2340 2341 2342 2343 2344 nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg %start_date)2345 2346 2347 2348 2349 2352 SUBROUTINE setup_variable_tables(ic_mode) 2353 CHARACTER(LEN=*), INTENT(IN) :: ic_mode 2354 INTEGER :: n_invar = 0 !< number of variables in the input variable table 2355 INTEGER :: n_outvar = 0 !< number of variables in the output variable table 2356 TYPE(nc_var), POINTER :: var 2357 2358 IF (TRIM(cfg%start_date) == '') THEN 2359 message = 'Simulation start date has not been set.' 2360 CALL inifor_abort('setup_variable_tables', message) 2361 ENDIF 2362 2363 nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg%start_date) 2364 2365 n_invar = 17 2366 n_outvar = 64 2367 ALLOCATE( input_var_table(n_invar) ) 2368 ALLOCATE( output_var_table(n_outvar) ) 2350 2369 2351 2370 ! … … 2353 2372 !- Section 1: NetCDF input variables 2354 2373 !------------------------------------------------------------------------------ 2355 2356 var %name = 'T_SO'2357 var %to_be_processed = .TRUE.2358 var %is_upside_down = .FALSE.2359 2360 2361 var %name = 'W_SO'2362 var %to_be_processed = .TRUE.2363 var %is_upside_down = .FALSE.2364 2365 2366 var %name = 'T'2367 var %to_be_processed = .TRUE.2368 var %is_upside_down = .TRUE.2369 2370 2371 var %name = 'QV'2372 var %to_be_processed = .TRUE.2373 var %is_upside_down = .TRUE.2374 2375 2376 var %name = 'U'2377 var %to_be_processed = .TRUE.2378 var %is_upside_down = .TRUE.2379 2380 2381 var %name = 'V'2382 var %to_be_processed = .TRUE.2383 var %is_upside_down = .TRUE.2384 2385 2386 var %name = 'W'2387 var %to_be_processed = .TRUE.2388 var %is_upside_down = .TRUE.2389 2390 2391 var %name = 'RAIN_GSP'2392 var %to_be_processed = .TRUE.2393 var %is_upside_down = .FALSE.2394 2395 2396 var %name = 'SNOW_GSP'2397 var %to_be_processed = .TRUE.2398 var %is_upside_down = .FALSE.2399 2400 2401 var %name = 'GRAU_GSP'2402 var %to_be_processed = .TRUE.2403 var %is_upside_down = .FALSE.2404 2405 2406 var %name = 'AEVAP_S'2407 var %to_be_processed = .TRUE.2408 var %is_upside_down = .FALSE.2409 2410 2411 var %name = 'T_2M'2412 var %to_be_processed = .TRUE.2413 var %is_upside_down = .FALSE.2414 2415 2416 var %name = 'ASWDIFD_S'2417 var %to_be_processed = .TRUE.2418 var %is_upside_down = .FALSE.2419 2420 2421 var %name = 'ASWDIR_S'2422 var %to_be_processed = .TRUE.2423 var %is_upside_down = .FALSE.2424 2425 2426 var %name = 'ASOB_S'2427 var %to_be_processed = .TRUE.2428 var %is_upside_down = .FALSE.2429 2430 2431 var %name = 'ATHB_S'2432 var %to_be_processed = .TRUE.2433 var %is_upside_down = .FALSE.2434 2435 2436 var %name = 'P'2437 var %to_be_processed = .TRUE.2438 var %is_upside_down = .TRUE.2374 var => input_var_table(1) 2375 var%name = 'T_SO' 2376 var%to_be_processed = .TRUE. 2377 var%is_upside_down = .FALSE. 2378 2379 var => input_var_table(2) 2380 var%name = 'W_SO' 2381 var%to_be_processed = .TRUE. 2382 var%is_upside_down = .FALSE. 2383 2384 var => input_var_table(3) 2385 var%name = 'T' 2386 var%to_be_processed = .TRUE. 2387 var%is_upside_down = .TRUE. 2388 2389 var => input_var_table(4) 2390 var%name = 'QV' 2391 var%to_be_processed = .TRUE. 2392 var%is_upside_down = .TRUE. 2393 2394 var => input_var_table(5) 2395 var%name = 'U' 2396 var%to_be_processed = .TRUE. 2397 var%is_upside_down = .TRUE. 2398 2399 var => input_var_table(6) 2400 var%name = 'V' 2401 var%to_be_processed = .TRUE. 2402 var%is_upside_down = .TRUE. 2403 2404 var => input_var_table(7) 2405 var%name = 'W' 2406 var%to_be_processed = .TRUE. 2407 var%is_upside_down = .TRUE. 2408 2409 var => input_var_table(8) 2410 var%name = 'RAIN_GSP' 2411 var%to_be_processed = .TRUE. 2412 var%is_upside_down = .FALSE. 2413 2414 var => input_var_table(9) 2415 var%name = 'SNOW_GSP' 2416 var%to_be_processed = .TRUE. 2417 var%is_upside_down = .FALSE. 2418 2419 var => input_var_table(10) 2420 var%name = 'GRAU_GSP' 2421 var%to_be_processed = .TRUE. 2422 var%is_upside_down = .FALSE. 2423 2424 var => input_var_table(11) 2425 var%name = 'AEVAP_S' 2426 var%to_be_processed = .TRUE. 2427 var%is_upside_down = .FALSE. 2428 2429 var => input_var_table(12) 2430 var%name = 'T_2M' 2431 var%to_be_processed = .TRUE. 2432 var%is_upside_down = .FALSE. 2433 2434 var => input_var_table(13) 2435 var%name = 'ASWDIFD_S' 2436 var%to_be_processed = .TRUE. 2437 var%is_upside_down = .FALSE. 2438 2439 var => input_var_table(14) 2440 var%name = 'ASWDIR_S' 2441 var%to_be_processed = .TRUE. 2442 var%is_upside_down = .FALSE. 2443 2444 var => input_var_table(15) 2445 var%name = 'ASOB_S' 2446 var%to_be_processed = .TRUE. 2447 var%is_upside_down = .FALSE. 2448 2449 var => input_var_table(16) 2450 var%name = 'ATHB_S' 2451 var%to_be_processed = .TRUE. 2452 var%is_upside_down = .FALSE. 2453 2454 var => input_var_table(17) 2455 var%name = 'P' 2456 var%to_be_processed = .TRUE. 2457 var%is_upside_down = .TRUE. 2439 2458 2440 2459 ! … … 2446 2465 ! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions 2447 2466 !------------------------------------------------------------------------------ 2448 2449 2450 2451 2452 2453 2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465 2466 2467 2468 2469 2470 2471 2472 2473 2474 2475 2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 output_var_table(3) %averaging_grid => averaged_initial_scalar_profile2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 2507 2508 2509 2510 2511 2512 2513 2514 2515 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2526 2527 2528 2529 2530 2531 2532 2533 2534 2535 2536 2537 2538 2539 2540 2541 2542 2543 2544 2545 2546 2547 2548 2549 2550 2551 2552 2553 2554 2555 2556 2557 2558 2559 2560 2561 output_var_table(9) %averaging_grid => averaged_initial_scalar_profile2562 2563 2564 2565 2566 2567 2568 2569 2570 2571 2572 2573 2574 2575 2576 2577 2578 2579 2580 2581 2582 2583 2584 2585 2586 2587 2588 2589 2590 2591 2592 2593 2594 2595 2596 2597 2598 2599 2600 2601 2602 2603 2604 2605 2606 2607 2608 2609 2610 2611 2612 2613 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2625 2626 2627 2628 2629 2630 2631 2632 2633 2634 2635 2636 2637 output_var_table(15) %averaging_grid => averaged_initial_scalar_profile2638 2639 2640 2641 2642 2643 2644 2645 2646 2647 2648 2649 2650 2651 2652 2653 2654 2655 2656 2657 2658 2659 2660 2661 2662 2663 2664 2665 2666 2667 2668 2669 2670 2671 2672 2673 2674 2675 2676 2677 2678 2679 2680 2681 2682 2683 2684 2685 2686 2687 2688 2689 2690 2691 2692 2693 2694 2695 2696 2697 2698 2699 2700 2701 2702 2703 2704 2705 2706 2707 2708 2709 2710 2711 2712 2713 output_var_table(21) %averaging_grid => averaged_initial_scalar_profile2714 2715 2716 2717 2718 2719 2720 2721 2722 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734 2735 2736 2737 2738 2739 2740 2741 2742 2743 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2759 2760 2761 2762 2763 2764 2765 2766 2767 2768 2769 2770 2771 2772 2773 2774 2775 2776 2777 2778 2779 2780 2781 2782 2783 2784 2785 2786 2787 2788 2789 output_var_table(27) %averaging_grid => averaged_initial_w_profile2790 2791 2792 2793 2794 2795 2796 2797 2798 2799 2800 2801 2802 2803 2804 2805 2806 2807 2808 2809 2810 2811 2812 2813 2814 2815 2816 2817 2818 2819 2820 2821 2822 2823 2824 2825 2826 2827 2828 2829 2830 2831 2832 2833 2834 2835 2836 2837 2838 2839 2840 2841 2842 2843 2844 2845 2846 2847 2848 2849 2850 2851 2852 2853 2854 2855 2856 2857 2858 2859 2860 2861 2862 2863 2864 2865 2866 2867 2868 2869 2870 2871 2872 2873 2874 2875 2876 2877 2878 2879 2880 2881 2882 2883 2884 2885 2886 2887 2888 2889 2890 2891 2892 2893 2894 2895 2896 2897 2898 2899 2900 2901 2902 2903 2904 2905 2906 2907 2908 2909 2910 2911 2912 2913 2914 2915 2916 2917 2918 2919 2920 2921 2922 2923 2924 2925 2926 2927 2928 2929 2930 2931 2932 2933 2934 2935 2936 2937 2938 2939 2940 2941 2942 2943 2944 2945 2946 2947 2948 2949 2950 2951 2952 2953 2954 2955 2956 2957 2958 2467 output_var_table(1) = init_nc_var( & 2468 name = 'init_soil_t', & 2469 std_name = "", & 2470 long_name = "initial soil temperature", & 2471 units = "K", & 2472 kind = "init soil", & 2473 input_id = 1, & 2474 output_file = output_file, & 2475 grid = palm_grid, & 2476 intermediate_grid = palm_intermediate & 2477 ) 2478 2479 output_var_table(2) = init_nc_var( & 2480 name = 'init_soil_m', & 2481 std_name = "", & 2482 long_name = "initial soil moisture", & 2483 units = "m^3/m^3", & 2484 kind = "init soil", & 2485 input_id = 1, & 2486 output_file = output_file, & 2487 grid = palm_grid, & 2488 intermediate_grid = palm_intermediate & 2489 ) 2490 2491 output_var_table(3) = init_nc_var( & 2492 name = 'init_atmosphere_pt', & 2493 std_name = "", & 2494 long_name = "initial potential temperature", & 2495 units = "K", & 2496 kind = "init scalar", & 2497 input_id = 1, & ! first in (T, p) IO group 2498 output_file = output_file, & 2499 grid = palm_grid, & 2500 intermediate_grid = palm_intermediate, & 2501 is_profile = (TRIM(ic_mode) == 'profile') & 2502 ) 2503 IF (TRIM(ic_mode) == 'profile') THEN 2504 output_var_table(3)%averaging_grid => averaged_initial_scalar_profile 2505 ENDIF 2506 2507 output_var_table(4) = init_nc_var( & 2508 name = 'ls_forcing_left_pt', & 2509 std_name = "", & 2510 long_name = "large-scale forcing for left model boundary for the potential temperature", & 2511 units = "K", & 2512 kind = "left scalar", & 2513 input_id = 1, & 2514 grid = scalars_west_grid, & 2515 intermediate_grid = scalars_west_intermediate, & 2516 output_file = output_file & 2517 ) 2518 2519 output_var_table(5) = init_nc_var( & 2520 name = 'ls_forcing_right_pt', & 2521 std_name = "", & 2522 long_name = "large-scale forcing for right model boundary for the potential temperature", & 2523 units = "K", & 2524 kind = "right scalar", & 2525 input_id = 1, & 2526 grid = scalars_east_grid, & 2527 intermediate_grid = scalars_east_intermediate, & 2528 output_file = output_file & 2529 ) 2530 2531 output_var_table(6) = init_nc_var( & 2532 name = 'ls_forcing_north_pt', & 2533 std_name = "", & 2534 long_name = "large-scale forcing for north model boundary for the potential temperature", & 2535 units = "K", & 2536 kind = "north scalar", & 2537 input_id = 1, & 2538 grid = scalars_north_grid, & 2539 intermediate_grid = scalars_north_intermediate, & 2540 output_file = output_file & 2541 ) 2542 2543 output_var_table(7) = init_nc_var( & 2544 name = 'ls_forcing_south_pt', & 2545 std_name = "", & 2546 long_name = "large-scale forcing for south model boundary for the potential temperature", & 2547 units = "K", & 2548 kind = "south scalar", & 2549 input_id = 1, & 2550 grid = scalars_south_grid, & 2551 intermediate_grid = scalars_south_intermediate, & 2552 output_file = output_file & 2553 ) 2554 2555 output_var_table(8) = init_nc_var( & 2556 name = 'ls_forcing_top_pt', & 2557 std_name = "", & 2558 long_name = "large-scale forcing for top model boundary for the potential temperature", & 2559 units = "K", & 2560 kind = "top scalar", & 2561 input_id = 1, & 2562 grid = scalars_top_grid, & 2563 intermediate_grid = scalars_top_intermediate, & 2564 output_file = output_file & 2565 ) 2566 2567 output_var_table(9) = init_nc_var( & 2568 name = 'init_atmosphere_qv', & 2569 std_name = "", & 2570 long_name = "initial specific humidity", & 2571 units = "kg/kg", & 2572 kind = "init scalar", & 2573 input_id = 3, & 2574 output_file = output_file, & 2575 grid = palm_grid, & 2576 intermediate_grid = palm_intermediate, & 2577 is_profile = (TRIM(ic_mode) == 'profile') & 2578 ) 2579 IF (TRIM(ic_mode) == 'profile') THEN 2580 output_var_table(9)%averaging_grid => averaged_initial_scalar_profile 2581 ENDIF 2582 2583 output_var_table(10) = init_nc_var( & 2584 name = 'ls_forcing_left_qv', & 2585 std_name = "", & 2586 long_name = "large-scale forcing for left model boundary for the specific humidity", & 2587 units = "kg/kg", & 2588 kind = "left scalar", & 2589 input_id = 3, & 2590 output_file = output_file, & 2591 grid = scalars_west_grid, & 2592 intermediate_grid = scalars_west_intermediate & 2593 ) 2594 2595 output_var_table(11) = init_nc_var( & 2596 name = 'ls_forcing_right_qv', & 2597 std_name = "", & 2598 long_name = "large-scale forcing for right model boundary for the specific humidity", & 2599 units = "kg/kg", & 2600 kind = "right scalar", & 2601 input_id = 3, & 2602 output_file = output_file, & 2603 grid = scalars_east_grid, & 2604 intermediate_grid = scalars_east_intermediate & 2605 ) 2606 2607 output_var_table(12) = init_nc_var( & 2608 name = 'ls_forcing_north_qv', & 2609 std_name = "", & 2610 long_name = "large-scale forcing for north model boundary for the specific humidity", & 2611 units = "kg/kg", & 2612 kind = "north scalar", & 2613 input_id = 3, & 2614 output_file = output_file, & 2615 grid = scalars_north_grid, & 2616 intermediate_grid = scalars_north_intermediate & 2617 ) 2618 2619 output_var_table(13) = init_nc_var( & 2620 name = 'ls_forcing_south_qv', & 2621 std_name = "", & 2622 long_name = "large-scale forcing for south model boundary for the specific humidity", & 2623 units = "kg/kg", & 2624 kind = "south scalar", & 2625 input_id = 3, & 2626 output_file = output_file, & 2627 grid = scalars_south_grid, & 2628 intermediate_grid = scalars_south_intermediate & 2629 ) 2630 2631 output_var_table(14) = init_nc_var( & 2632 name = 'ls_forcing_top_qv', & 2633 std_name = "", & 2634 long_name = "large-scale forcing for top model boundary for the specific humidity", & 2635 units = "kg/kg", & 2636 kind = "top scalar", & 2637 input_id = 3, & 2638 output_file = output_file, & 2639 grid = scalars_top_grid, & 2640 intermediate_grid = scalars_top_intermediate & 2641 ) 2642 2643 output_var_table(15) = init_nc_var( & 2644 name = 'init_atmosphere_u', & 2645 std_name = "", & 2646 long_name = "initial wind component in x direction", & 2647 units = "m/s", & 2648 kind = "init u", & 2649 input_id = 1, & ! first in (U, V) I/O group 2650 output_file = output_file, & 2651 grid = u_initial_grid, & 2652 intermediate_grid = u_initial_intermediate, & 2653 is_profile = (TRIM(ic_mode) == 'profile') & 2654 ) 2655 IF (TRIM(ic_mode) == 'profile') THEN 2656 output_var_table(15)%averaging_grid => averaged_initial_scalar_profile 2657 ENDIF 2658 2659 output_var_table(16) = init_nc_var( & 2660 name = 'ls_forcing_left_u', & 2661 std_name = "", & 2662 long_name = "large-scale forcing for left model boundary for the wind component in x direction", & 2663 units = "m/s", & 2664 kind = "left u", & 2665 input_id = 1, & ! first in (U, V) I/O group 2666 output_file = output_file, & 2667 grid = u_west_grid, & 2668 intermediate_grid = u_west_intermediate & 2669 ) 2670 2671 output_var_table(17) = init_nc_var( & 2672 name = 'ls_forcing_right_u', & 2673 std_name = "", & 2674 long_name = "large-scale forcing for right model boundary for the wind component in x direction", & 2675 units = "m/s", & 2676 kind = "right u", & 2677 input_id = 1, & ! first in (U, V) I/O group 2678 output_file = output_file, & 2679 grid = u_east_grid, & 2680 intermediate_grid = u_east_intermediate & 2681 ) 2682 2683 output_var_table(18) = init_nc_var( & 2684 name = 'ls_forcing_north_u', & 2685 std_name = "", & 2686 long_name = "large-scale forcing for north model boundary for the wind component in x direction", & 2687 units = "m/s", & 2688 kind = "north u", & 2689 input_id = 1, & ! first in (U, V) I/O group 2690 output_file = output_file, & 2691 grid = u_north_grid, & 2692 intermediate_grid = u_north_intermediate & 2693 ) 2694 2695 output_var_table(19) = init_nc_var( & 2696 name = 'ls_forcing_south_u', & 2697 std_name = "", & 2698 long_name = "large-scale forcing for south model boundary for the wind component in x direction", & 2699 units = "m/s", & 2700 kind = "south u", & 2701 input_id = 1, & ! first in (U, V) I/O group 2702 output_file = output_file, & 2703 grid = u_south_grid, & 2704 intermediate_grid = u_south_intermediate & 2705 ) 2706 2707 output_var_table(20) = init_nc_var( & 2708 name = 'ls_forcing_top_u', & 2709 std_name = "", & 2710 long_name = "large-scale forcing for top model boundary for the wind component in x direction", & 2711 units = "m/s", & 2712 kind = "top u", & 2713 input_id = 1, & ! first in (U, V) I/O group 2714 output_file = output_file, & 2715 grid = u_top_grid, & 2716 intermediate_grid = u_top_intermediate & 2717 ) 2718 2719 output_var_table(21) = init_nc_var( & 2720 name = 'init_atmosphere_v', & 2721 std_name = "", & 2722 long_name = "initial wind component in y direction", & 2723 units = "m/s", & 2724 kind = "init v", & 2725 input_id = 2, & ! second in (U, V) I/O group 2726 output_file = output_file, & 2727 grid = v_initial_grid, & 2728 intermediate_grid = v_initial_intermediate, & 2729 is_profile = (TRIM(ic_mode) == 'profile') & 2730 ) 2731 IF (TRIM(ic_mode) == 'profile') THEN 2732 output_var_table(21)%averaging_grid => averaged_initial_scalar_profile 2733 ENDIF 2734 2735 output_var_table(22) = init_nc_var( & 2736 name = 'ls_forcing_left_v', & 2737 std_name = "", & 2738 long_name = "large-scale forcing for left model boundary for the wind component in y direction", & 2739 units = "m/s", & 2740 kind = "right v", & 2741 input_id = 2, & ! second in (U, V) I/O group 2742 output_file = output_file, & 2743 grid = v_west_grid, & 2744 intermediate_grid = v_west_intermediate & 2745 ) 2746 2747 output_var_table(23) = init_nc_var( & 2748 name = 'ls_forcing_right_v', & 2749 std_name = "", & 2750 long_name = "large-scale forcing for right model boundary for the wind component in y direction", & 2751 units = "m/s", & 2752 kind = "right v", & 2753 input_id = 2, & ! second in (U, V) I/O group 2754 output_file = output_file, & 2755 grid = v_east_grid, & 2756 intermediate_grid = v_east_intermediate & 2757 ) 2758 2759 output_var_table(24) = init_nc_var( & 2760 name = 'ls_forcing_north_v', & 2761 std_name = "", & 2762 long_name = "large-scale forcing for north model boundary for the wind component in y direction", & 2763 units = "m/s", & 2764 kind = "north v", & 2765 input_id = 2, & ! second in (U, V) I/O group 2766 output_file = output_file, & 2767 grid = v_north_grid, & 2768 intermediate_grid = v_north_intermediate & 2769 ) 2770 2771 output_var_table(25) = init_nc_var( & 2772 name = 'ls_forcing_south_v', & 2773 std_name = "", & 2774 long_name = "large-scale forcing for south model boundary for the wind component in y direction", & 2775 units = "m/s", & 2776 kind = "south v", & 2777 input_id = 2, & ! second in (U, V) I/O group 2778 output_file = output_file, & 2779 grid = v_south_grid, & 2780 intermediate_grid = v_south_intermediate & 2781 ) 2782 2783 output_var_table(26) = init_nc_var( & 2784 name = 'ls_forcing_top_v', & 2785 std_name = "", & 2786 long_name = "large-scale forcing for top model boundary for the wind component in y direction", & 2787 units = "m/s", & 2788 kind = "top v", & 2789 input_id = 2, & ! second in (U, V) I/O group 2790 output_file = output_file, & 2791 grid = v_top_grid, & 2792 intermediate_grid = v_top_intermediate & 2793 ) 2794 2795 output_var_table(27) = init_nc_var( & 2796 name = 'init_atmosphere_w', & 2797 std_name = "", & 2798 long_name = "initial wind component in z direction", & 2799 units = "m/s", & 2800 kind = "init w", & 2801 input_id = 1, & 2802 output_file = output_file, & 2803 grid = w_initial_grid, & 2804 intermediate_grid = w_initial_intermediate, & 2805 is_profile = (TRIM(ic_mode) == 'profile') & 2806 ) 2807 IF (TRIM(ic_mode) == 'profile') THEN 2808 output_var_table(27)%averaging_grid => averaged_initial_w_profile 2809 ENDIF 2810 2811 output_var_table(28) = init_nc_var( & 2812 name = 'ls_forcing_left_w', & 2813 std_name = "", & 2814 long_name = "large-scale forcing for left model boundary for the wind component in z direction", & 2815 units = "m/s", & 2816 kind = "left w", & 2817 input_id = 1, & 2818 output_file = output_file, & 2819 grid = w_west_grid, & 2820 intermediate_grid = w_west_intermediate & 2821 ) 2822 2823 output_var_table(29) = init_nc_var( & 2824 name = 'ls_forcing_right_w', & 2825 std_name = "", & 2826 long_name = "large-scale forcing for right model boundary for the wind component in z direction", & 2827 units = "m/s", & 2828 kind = "right w", & 2829 input_id = 1, & 2830 output_file = output_file, & 2831 grid = w_east_grid, & 2832 intermediate_grid = w_east_intermediate & 2833 ) 2834 2835 output_var_table(30) = init_nc_var( & 2836 name = 'ls_forcing_north_w', & 2837 std_name = "", & 2838 long_name = "large-scale forcing for north model boundary for the wind component in z direction", & 2839 units = "m/s", & 2840 kind = "north w", & 2841 input_id = 1, & 2842 output_file = output_file, & 2843 grid = w_north_grid, & 2844 intermediate_grid = w_north_intermediate & 2845 ) 2846 2847 output_var_table(31) = init_nc_var( & 2848 name = 'ls_forcing_south_w', & 2849 std_name = "", & 2850 long_name = "large-scale forcing for south model boundary for the wind component in z direction", & 2851 units = "m/s", & 2852 kind = "south w", & 2853 input_id = 1, & 2854 output_file = output_file, & 2855 grid = w_south_grid, & 2856 intermediate_grid = w_south_intermediate & 2857 ) 2858 2859 output_var_table(32) = init_nc_var( & 2860 name = 'ls_forcing_top_w', & 2861 std_name = "", & 2862 long_name = "large-scale forcing for top model boundary for the wind component in z direction", & 2863 units = "m/s", & 2864 kind = "top w", & 2865 input_id = 1, & 2866 output_file = output_file, & 2867 grid = w_top_grid, & 2868 intermediate_grid = w_top_intermediate & 2869 ) 2870 2871 output_var_table(33) = init_nc_var( & 2872 name = 'ls_forcing_soil_rain', & 2873 std_name = "", & 2874 long_name = "large-scale forcing rain", & 2875 units = "kg/m2", & 2876 kind = "surface forcing", & 2877 input_id = 1, & 2878 output_file = output_file, & 2879 grid = palm_grid, & 2880 intermediate_grid = palm_intermediate & 2881 ) 2882 2883 output_var_table(34) = init_nc_var( & 2884 name = 'ls_forcing_soil_snow', & 2885 std_name = "", & 2886 long_name = "large-scale forcing snow", & 2887 units = "kg/m2", & 2888 kind = "surface forcing", & 2889 input_id = 1, & 2890 output_file = output_file, & 2891 grid = palm_grid, & 2892 intermediate_grid = palm_intermediate & 2893 ) 2894 2895 output_var_table(35) = init_nc_var( & 2896 name = 'ls_forcing_soil_graupel', & 2897 std_name = "", & 2898 long_name = "large-scale forcing graupel", & 2899 units = "kg/m2", & 2900 kind = "surface forcing", & 2901 input_id = 1, & 2902 output_file = output_file, & 2903 grid = palm_grid, & 2904 intermediate_grid = palm_intermediate & 2905 ) 2906 2907 output_var_table(36) = init_nc_var( & 2908 name = 'ls_forcing_soil_t_2m', & 2909 std_name = "", & 2910 long_name = "large-scale forcing 2m air temperature", & 2911 units = "kg/m2", & 2912 kind = "surface forcing", & 2913 input_id = 1, & 2914 output_file = output_file, & 2915 grid = palm_grid, & 2916 intermediate_grid = palm_intermediate & 2917 ) 2918 2919 output_var_table(37) = init_nc_var( & 2920 name = 'ls_forcing_soil_evap', & 2921 std_name = "", & 2922 long_name = "large-scale forcing evapo-transpiration", & 2923 units = "kg/m2", & 2924 kind = "surface forcing", & 2925 input_id = 1, & 2926 output_file = output_file, & 2927 grid = palm_grid, & 2928 intermediate_grid = palm_intermediate & 2929 ) 2930 2931 output_var_table(38) = init_nc_var( & 2932 name = 'rad_swd_dif_0', & 2933 std_name = "", & 2934 long_name = "incoming diffuse shortwave radiative flux at the surface", & 2935 units = "W/m2", & 2936 kind = "surface forcing", & 2937 input_id = 1, & 2938 output_file = output_file, & 2939 grid = palm_grid, & 2940 intermediate_grid = palm_intermediate & 2941 ) 2942 2943 output_var_table(39) = init_nc_var( & 2944 name = 'rad_swd_dir_0', & 2945 std_name = "", & 2946 long_name = "incoming direct shortwave radiative flux at the surface", & 2947 units = "W/m2", & 2948 kind = "surface forcing", & 2949 input_id = 1, & 2950 output_file = output_file, & 2951 grid = palm_grid, & 2952 intermediate_grid = palm_intermediate & 2953 ) 2954 2955 output_var_table(40) = init_nc_var( & 2956 name = 'rad_sw_bal_0', & 2957 std_name = "", & 2958 long_name = "shortwave radiation balance at the surface", & 2959 units = "W/m2", & 2960 kind = "surface forcing", & 2961 input_id = 1, & 2962 output_file = output_file, & 2963 grid = palm_grid, & 2964 intermediate_grid = palm_intermediate & 2965 ) 2966 2967 output_var_table(41) = init_nc_var( & 2968 name = 'rad_lw_bal_0', & 2969 std_name = "", & 2970 long_name = "longwave radiation balance at the surface", & 2971 units = "W/m2", & 2972 kind = "surface forcing", & 2973 input_id = 1, & 2974 output_file = output_file, & 2975 grid = palm_grid, & 2976 intermediate_grid = palm_intermediate & 2977 ) 2959 2978 ! 2960 2979 !------------------------------------------------------------------------------ 2961 2980 ! Section 2.2: Idealized large-scale forcings 2962 2981 !------------------------------------------------------------------------------ 2963 2964 2965 2966 2967 2968 2969 2970 2971 2972 2973 2974 output_var_table(42) %averaging_grid => averaged_scalar_profile2975 2976 2977 2978 2979 2980 2981 2982 2983 2984 2985 2986 2987 2988 2989 2990 2991 2992 2993 2994 2995 2996 2997 2998 2999 3000 3001 3002 3003 3004 3005 3006 3007 3008 3009 3010 3011 output_var_table(45) %to_be_processed = ls_forcing_variables_required3012 3013 3014 3015 3016 3017 3018 3019 3020 3021 3022 3023 3024 output_var_table(46) %to_be_processed = ls_forcing_variables_required3025 3026 3027 3028 3029 3030 3031 3032 3033 3034 3035 3036 3037 output_var_table(47) %to_be_processed = ls_forcing_variables_required3038 3039 3040 3041 3042 3043 3044 3045 3046 3047 3048 3049 3050 output_var_table(48) %to_be_processed = ls_forcing_variables_required3051 3052 3053 3054 3055 3056 3057 3058 3059 3060 3061 3062 3063 3064 output_var_table(49) %to_be_processed = ls_forcing_variables_required3065 3066 3067 3068 3069 3070 3071 3072 3073 3074 3075 3076 3077 output_var_table(50) %to_be_processed = ls_forcing_variables_required3078 3079 3080 3081 3082 3083 3084 3085 3086 3087 3088 3089 3090 output_var_table(51) %to_be_processed = ls_forcing_variables_required3091 3092 3093 3094 3095 3096 3097 3098 3099 3100 3101 3102 3103 output_var_table(52) %to_be_processed = ls_forcing_variables_required3104 3105 3106 3107 3108 3109 3110 3111 3112 3113 3114 3115 3116 3117 output_var_table(53) %to_be_processed = ls_forcing_variables_required3118 3119 3120 3121 3122 3123 3124 3125 3126 3127 3128 3129 3130 output_var_table(54) %to_be_processed = ls_forcing_variables_required3131 3132 3133 3134 3135 3136 3137 3138 3139 3140 3141 3142 3143 output_var_table(55) %to_be_processed = ls_forcing_variables_required3144 3145 3146 3147 3148 3149 3150 3151 3152 3153 3154 3155 3156 3157 output_var_table(56) %averaging_grid => averaged_scalar_profile3158 3159 3160 3161 3162 3163 3164 3165 3166 3167 3168 3169 3170 3171 output_var_table(57) %averaging_grid => north_averaged_scalar_profile3172 output_var_table(57) % to_be_processed = .NOT. cfg %ug_defined_by_user3173 3174 3175 3176 3177 3178 3179 3180 3181 3182 3183 3184 3185 3186 output_var_table(58) %averaging_grid => south_averaged_scalar_profile3187 output_var_table(58) % to_be_processed = .NOT. cfg %ug_defined_by_user3188 3189 3190 3191 3192 3193 3194 3195 3196 3197 3198 3199 3200 3201 output_var_table(59) %averaging_grid => east_averaged_scalar_profile3202 output_var_table(59) % to_be_processed = .NOT. cfg %ug_defined_by_user3203 3204 3205 3206 3207 3208 3209 3210 3211 3212 3213 3214 3215 3216 output_var_table(60) %averaging_grid => west_averaged_scalar_profile3217 output_var_table(60) % to_be_processed = .NOT. cfg %ug_defined_by_user3218 3219 3220 3221 3222 3223 3224 3225 3226 3227 3228 3229 3230 output_var_table(61) %averaging_grid => north_averaged_scalar_profile3231 output_var_table(61) % to_be_processed = .NOT. cfg %ug_defined_by_user3232 3233 3234 3235 3236 3237 3238 3239 3240 3241 3242 3243 3244 3245 output_var_table(62) %averaging_grid => south_averaged_scalar_profile3246 output_var_table(62) % to_be_processed = .NOT. cfg %ug_defined_by_user3247 3248 3249 3250 3251 3252 3253 3254 3255 3256 3257 3258 3259 3260 output_var_table(63) %averaging_grid => east_averaged_scalar_profile3261 output_var_table(63) % to_be_processed = .NOT. cfg %ug_defined_by_user3262 3263 3264 3265 3266 3267 3268 3269 3270 3271 3272 3273 3274 3275 output_var_table(64) %averaging_grid => west_averaged_scalar_profile3276 output_var_table(64) % to_be_processed = .NOT. cfg %ug_defined_by_user3277 3278 ! 3279 !-- 3280 output_var_table(:) %source = nc_source_text3281 3282 3283 2982 output_var_table(42) = init_nc_var( & 2983 name = 'surface_forcing_surface_pressure', & 2984 std_name = "", & 2985 long_name = "surface pressure", & 2986 units = "Pa", & 2987 kind = "time series", & 2988 input_id = 2, & ! second in (T, p) I/O group 2989 output_file = output_file, & 2990 grid = palm_grid, & 2991 intermediate_grid = palm_intermediate & 2992 ) 2993 output_var_table(42)%averaging_grid => averaged_scalar_profile 2994 2995 output_var_table(43) = init_nc_var( & 2996 name = 'ls_forcing_ug', & 2997 std_name = "", & 2998 long_name = "geostrophic wind (u component)", & 2999 units = "m/s", & 3000 kind = "geostrophic", & 3001 input_id = 1, & 3002 output_file = output_file, & 3003 grid = averaged_scalar_profile, & 3004 intermediate_grid = averaged_scalar_profile & 3005 ) 3006 3007 output_var_table(44) = init_nc_var( & 3008 name = 'ls_forcing_vg', & 3009 std_name = "", & 3010 long_name = "geostrophic wind (v component)", & 3011 units = "m/s", & 3012 kind = "geostrophic", & 3013 input_id = 1, & 3014 output_file = output_file, & 3015 grid = averaged_scalar_profile, & 3016 intermediate_grid = averaged_scalar_profile & 3017 ) 3018 3019 output_var_table(45) = init_nc_var( & 3020 name = 'nudging_u', & 3021 std_name = "", & 3022 long_name = "wind component in x direction", & 3023 units = "m/s", & 3024 kind = "geostrophic", & 3025 input_id = 1, & 3026 output_file = output_file, & 3027 grid = averaged_scalar_profile, & 3028 intermediate_grid = averaged_scalar_profile & 3029 ) 3030 output_var_table(45)%to_be_processed = ls_forcing_variables_required 3031 3032 output_var_table(46) = init_nc_var( & 3033 name = 'nudging_v', & 3034 std_name = "", & 3035 long_name = "wind component in y direction", & 3036 units = "m/s", & 3037 kind = "large-scale scalar forcing", & 3038 input_id = 1, & 3039 output_file = output_file, & 3040 grid = averaged_scalar_profile, & 3041 intermediate_grid = averaged_scalar_profile & 3042 ) 3043 output_var_table(46)%to_be_processed = ls_forcing_variables_required 3044 3045 output_var_table(47) = init_nc_var( & 3046 name = 'ls_forcing_sub_w', & 3047 std_name = "", & 3048 long_name = "subsidence velocity of w", & 3049 units = "m/s", & 3050 kind = "large-scale w forcing", & 3051 input_id = 1, & 3052 output_file = output_file, & 3053 grid = averaged_scalar_profile, & 3054 intermediate_grid = averaged_scalar_profile & 3055 ) 3056 output_var_table(47)%to_be_processed = ls_forcing_variables_required 3057 3058 output_var_table(48) = init_nc_var( & 3059 name = 'nudging_w', & 3060 std_name = "", & 3061 long_name = "wind component in w direction", & 3062 units = "m/s", & 3063 kind = "large-scale w forcing", & 3064 input_id = 1, & 3065 output_file = output_file, & 3066 grid = averaged_w_profile, & 3067 intermediate_grid = averaged_w_profile & 3068 ) 3069 output_var_table(48)%to_be_processed = ls_forcing_variables_required 3070 3071 3072 output_var_table(49) = init_nc_var( & 3073 name = 'ls_forcing_adv_pt', & 3074 std_name = "", & 3075 long_name = "advection of potential temperature", & 3076 units = "K/s", & 3077 kind = "large-scale scalar forcing", & 3078 input_id = 1, & 3079 output_file = output_file, & 3080 grid = averaged_scalar_profile, & 3081 intermediate_grid = averaged_scalar_profile & 3082 ) 3083 output_var_table(49)%to_be_processed = ls_forcing_variables_required 3084 3085 output_var_table(50) = init_nc_var( & 3086 name = 'ls_forcing_sub_pt', & 3087 std_name = "", & 3088 long_name = "subsidence velocity of potential temperature", & 3089 units = "K/s", & 3090 kind = "large-scale scalar forcing", & 3091 input_id = 1, & 3092 output_file = output_file, & 3093 grid = averaged_scalar_profile, & 3094 intermediate_grid = averaged_scalar_profile & 3095 ) 3096 output_var_table(50)%to_be_processed = ls_forcing_variables_required 3097 3098 output_var_table(51) = init_nc_var( & 3099 name = 'nudging_pt', & 3100 std_name = "", & 3101 long_name = "potential temperature", & 3102 units = "K", & 3103 kind = "large-scale scalar forcing", & 3104 input_id = 1, & 3105 output_file = output_file, & 3106 grid = averaged_scalar_profile, & 3107 intermediate_grid = averaged_scalar_profile & 3108 ) 3109 output_var_table(51)%to_be_processed = ls_forcing_variables_required 3110 3111 output_var_table(52) = init_nc_var( & 3112 name = 'ls_forcing_adv_qv', & 3113 std_name = "", & 3114 long_name = "advection of specific humidity", & 3115 units = "kg/kg/s", & 3116 kind = "large-scale scalar forcing", & 3117 input_id = 3, & 3118 output_file = output_file, & 3119 grid = averaged_scalar_profile, & 3120 intermediate_grid = averaged_scalar_profile & 3121 ) 3122 output_var_table(52)%to_be_processed = ls_forcing_variables_required 3123 3124 3125 output_var_table(53) = init_nc_var( & 3126 name = 'ls_forcing_sub_qv', & 3127 std_name = "", & 3128 long_name = "subsidence velocity of specific humidity", & 3129 units = "kg/kg/s", & 3130 kind = "large-scale scalar forcing", & 3131 input_id = 3, & 3132 output_file = output_file, & 3133 grid = averaged_scalar_profile, & 3134 intermediate_grid = averaged_scalar_profile & 3135 ) 3136 output_var_table(53)%to_be_processed = ls_forcing_variables_required 3137 3138 output_var_table(54) = init_nc_var( & 3139 name = 'nudging_qv', & 3140 std_name = "", & 3141 long_name = "specific humidity", & 3142 units = "kg/kg", & 3143 kind = "large-scale scalar forcing", & 3144 input_id = 3, & 3145 output_file = output_file, & 3146 grid = averaged_scalar_profile, & 3147 intermediate_grid = averaged_scalar_profile & 3148 ) 3149 output_var_table(54)%to_be_processed = ls_forcing_variables_required 3150 3151 output_var_table(55) = init_nc_var( & 3152 name = 'nudging_tau', & 3153 std_name = "", & 3154 long_name = "nudging relaxation time scale", & 3155 units = "s", & 3156 kind = "constant scalar profile", & 3157 input_id = 1, & 3158 output_file = output_file, & 3159 grid = averaged_scalar_profile, & 3160 intermediate_grid = averaged_scalar_profile & 3161 ) 3162 output_var_table(55)%to_be_processed = ls_forcing_variables_required 3163 3164 3165 output_var_table(56) = init_nc_var( & 3166 name = 'internal_density_centre', & 3167 std_name = "", & 3168 long_name = "", & 3169 units = "", & 3170 kind = "internal profile", & 3171 input_id = 4, & 3172 output_file = output_file, & 3173 grid = averaged_scalar_profile, & 3174 intermediate_grid = averaged_scalar_profile & 3175 ) 3176 output_var_table(56)%averaging_grid => averaged_scalar_profile 3177 3178 3179 output_var_table(57) = init_nc_var( & 3180 name = 'internal_density_north', & 3181 std_name = "", & 3182 long_name = "", & 3183 units = "", & 3184 kind = "internal profile", & 3185 input_id = 4, & 3186 output_file = output_file, & 3187 grid = north_averaged_scalar_profile, & 3188 intermediate_grid = north_averaged_scalar_profile & 3189 ) 3190 output_var_table(57)%averaging_grid => north_averaged_scalar_profile 3191 output_var_table(57)%to_be_processed = .NOT. cfg%ug_defined_by_user 3192 3193 3194 output_var_table(58) = init_nc_var( & 3195 name = 'internal_density_south', & 3196 std_name = "", & 3197 long_name = "", & 3198 units = "", & 3199 kind = "internal profile", & 3200 input_id = 4, & 3201 output_file = output_file, & 3202 grid = south_averaged_scalar_profile, & 3203 intermediate_grid = south_averaged_scalar_profile & 3204 ) 3205 output_var_table(58)%averaging_grid => south_averaged_scalar_profile 3206 output_var_table(58)%to_be_processed = .NOT. cfg%ug_defined_by_user 3207 3208 3209 output_var_table(59) = init_nc_var( & 3210 name = 'internal_density_east', & 3211 std_name = "", & 3212 long_name = "", & 3213 units = "", & 3214 kind = "internal profile", & 3215 input_id = 4, & 3216 output_file = output_file, & 3217 grid = east_averaged_scalar_profile, & 3218 intermediate_grid = east_averaged_scalar_profile & 3219 ) 3220 output_var_table(59)%averaging_grid => east_averaged_scalar_profile 3221 output_var_table(59)%to_be_processed = .NOT. cfg%ug_defined_by_user 3222 3223 3224 output_var_table(60) = init_nc_var( & 3225 name = 'internal_density_west', & 3226 std_name = "", & 3227 long_name = "", & 3228 units = "", & 3229 kind = "internal profile", & 3230 input_id = 4, & 3231 output_file = output_file, & 3232 grid = west_averaged_scalar_profile, & 3233 intermediate_grid = west_averaged_scalar_profile & 3234 ) 3235 output_var_table(60)%averaging_grid => west_averaged_scalar_profile 3236 output_var_table(60)%to_be_processed = .NOT. cfg%ug_defined_by_user 3237 3238 output_var_table(61) = init_nc_var( & 3239 name = 'internal_pressure_north', & 3240 std_name = "", & 3241 long_name = "", & 3242 units = "", & 3243 kind = "internal profile", & 3244 input_id = 2, & 3245 output_file = output_file, & 3246 grid = north_averaged_scalar_profile, & 3247 intermediate_grid = north_averaged_scalar_profile & 3248 ) 3249 output_var_table(61)%averaging_grid => north_averaged_scalar_profile 3250 output_var_table(61)%to_be_processed = .NOT. cfg%ug_defined_by_user 3251 3252 3253 output_var_table(62) = init_nc_var( & 3254 name = 'internal_pressure_south', & 3255 std_name = "", & 3256 long_name = "", & 3257 units = "", & 3258 kind = "internal profile", & 3259 input_id = 2, & 3260 output_file = output_file, & 3261 grid = south_averaged_scalar_profile, & 3262 intermediate_grid = south_averaged_scalar_profile & 3263 ) 3264 output_var_table(62)%averaging_grid => south_averaged_scalar_profile 3265 output_var_table(62)%to_be_processed = .NOT. cfg%ug_defined_by_user 3266 3267 3268 output_var_table(63) = init_nc_var( & 3269 name = 'internal_pressure_east', & 3270 std_name = "", & 3271 long_name = "", & 3272 units = "", & 3273 kind = "internal profile", & 3274 input_id = 2, & 3275 output_file = output_file, & 3276 grid = east_averaged_scalar_profile, & 3277 intermediate_grid = east_averaged_scalar_profile & 3278 ) 3279 output_var_table(63)%averaging_grid => east_averaged_scalar_profile 3280 output_var_table(63)%to_be_processed = .NOT. cfg%ug_defined_by_user 3281 3282 3283 output_var_table(64) = init_nc_var( & 3284 name = 'internal_pressure_west', & 3285 std_name = "", & 3286 long_name = "", & 3287 units = "", & 3288 kind = "internal profile", & 3289 input_id = 2, & 3290 output_file = output_file, & 3291 grid = west_averaged_scalar_profile, & 3292 intermediate_grid = west_averaged_scalar_profile & 3293 ) 3294 output_var_table(64)%averaging_grid => west_averaged_scalar_profile 3295 output_var_table(64)%to_be_processed = .NOT. cfg%ug_defined_by_user 3296 3297 ! 3298 !-- Attributes shared among all variables 3299 output_var_table(:)%source = nc_source_text 3300 3301 3302 END SUBROUTINE setup_variable_tables 3284 3303 3285 3304 … … 3291 3310 !> 'lod', as defined by the PALM-4U input data standard. 3292 3311 !------------------------------------------------------------------------------! 3293 3294 3295 3296 3297 3298 3299 3300 3301 3302 3303 3304 3305 3306 3307 3308 3309 3310 3311 3312 var %name = name3313 var %standard_name = std_name3314 var %long_name = long_name3315 var %units = units3316 var %kind = TRIM(out_var_kind)3317 var %input_id = input_id3318 var % nt = SIZE (output_file %time)3319 var %grid => grid3320 var %intermediate_grid => intermediate_grid3321 3322 3312 FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, & 3313 grid, intermediate_grid, output_file, is_profile) & 3314 RESULT(var) 3315 3316 CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind 3317 INTEGER, INTENT(IN) :: input_id 3318 TYPE(grid_definition), INTENT(IN), TARGET :: grid, intermediate_grid 3319 TYPE(nc_file), INTENT(IN) :: output_file 3320 LOGICAL, INTENT(IN), OPTIONAL :: is_profile 3321 3322 CHARACTER(LEN=LNAME) :: out_var_kind 3323 TYPE(nc_var) :: var 3324 3325 out_var_kind = TRIM(kind) 3326 3327 IF (PRESENT(is_profile)) THEN 3328 IF (is_profile) out_var_kind = TRIM(kind) // ' profile' 3329 ENDIF 3330 3331 var%name = name 3332 var%standard_name = std_name 3333 var%long_name = long_name 3334 var%units = units 3335 var%kind = TRIM(out_var_kind) 3336 var%input_id = input_id 3337 var%nt = SIZE (output_file%time) 3338 var%grid => grid 3339 var%intermediate_grid => intermediate_grid 3340 3341 SELECT CASE( TRIM(out_var_kind) ) 3323 3342 3324 3343 ! … … 3327 3346 !-- TODO: and pass into init_nc_var. 3328 3347 CASE( 'init soil' ) 3329 var %nt = 13330 var %lod = 23331 var %ndim = 33332 var % dimids(1:3) = output_file %dimids_soil3333 var % dimvarids(1:3) = output_file %dimvarids_soil3334 var %to_be_processed = init_variables_required3335 var %is_internal = .FALSE.3336 var %task = "interpolate_2d"3348 var%nt = 1 3349 var%lod = 2 3350 var%ndim = 3 3351 var%dimids(1:3) = output_file%dimids_soil 3352 var%dimvarids(1:3) = output_file%dimvarids_soil 3353 var%to_be_processed = init_variables_required 3354 var%is_internal = .FALSE. 3355 var%task = "interpolate_2d" 3337 3356 3338 3357 CASE( 'init scalar' ) 3339 var %nt = 13340 var %lod = 23341 var %ndim = 33342 var % dimids(1:3) = output_file %dimids_scl3343 var % dimvarids(1:3) = output_file %dimvarids_scl3344 var %to_be_processed = init_variables_required3345 var %is_internal = .FALSE.3346 var %task = "interpolate_3d"3358 var%nt = 1 3359 var%lod = 2 3360 var%ndim = 3 3361 var%dimids(1:3) = output_file%dimids_scl 3362 var%dimvarids(1:3) = output_file%dimvarids_scl 3363 var%to_be_processed = init_variables_required 3364 var%is_internal = .FALSE. 3365 var%task = "interpolate_3d" 3347 3366 3348 3367 CASE( 'init u' ) 3349 var %nt = 13350 var %lod = 23351 var %ndim = 33352 var % dimids(1) = output_file %dimids_vel(1)3353 var % dimids(2) = output_file %dimids_scl(2)3354 var % dimids(3) = output_file %dimids_scl(3)3355 var % dimvarids(1) = output_file %dimvarids_vel(1)3356 var % dimvarids(2) = output_file %dimvarids_scl(2)3357 var % dimvarids(3) = output_file %dimvarids_scl(3)3358 var %to_be_processed = init_variables_required3359 var %is_internal = .FALSE.3360 var %task = "interpolate_3d"3368 var%nt = 1 3369 var%lod = 2 3370 var%ndim = 3 3371 var%dimids(1) = output_file%dimids_vel(1) 3372 var%dimids(2) = output_file%dimids_scl(2) 3373 var%dimids(3) = output_file%dimids_scl(3) 3374 var%dimvarids(1) = output_file%dimvarids_vel(1) 3375 var%dimvarids(2) = output_file%dimvarids_scl(2) 3376 var%dimvarids(3) = output_file%dimvarids_scl(3) 3377 var%to_be_processed = init_variables_required 3378 var%is_internal = .FALSE. 3379 var%task = "interpolate_3d" 3361 3380 3362 3381 CASE( 'init v' ) 3363 var %nt = 13364 var %lod = 23365 var %ndim = 33366 var % dimids(1) = output_file %dimids_scl(1)3367 var % dimids(2) = output_file %dimids_vel(2)3368 var % dimids(3) = output_file %dimids_scl(3)3369 var % dimvarids(1) = output_file %dimvarids_scl(1)3370 var % dimvarids(2) = output_file %dimvarids_vel(2)3371 var % dimvarids(3) = output_file %dimvarids_scl(3)3372 var %to_be_processed = init_variables_required3373 var %is_internal = .FALSE.3374 var %task = "interpolate_3d"3382 var%nt = 1 3383 var%lod = 2 3384 var%ndim = 3 3385 var%dimids(1) = output_file%dimids_scl(1) 3386 var%dimids(2) = output_file%dimids_vel(2) 3387 var%dimids(3) = output_file%dimids_scl(3) 3388 var%dimvarids(1) = output_file%dimvarids_scl(1) 3389 var%dimvarids(2) = output_file%dimvarids_vel(2) 3390 var%dimvarids(3) = output_file%dimvarids_scl(3) 3391 var%to_be_processed = init_variables_required 3392 var%is_internal = .FALSE. 3393 var%task = "interpolate_3d" 3375 3394 3376 3395 CASE( 'init w' ) 3377 var %nt = 13378 var %lod = 23379 var %ndim = 33380 var % dimids(1) = output_file %dimids_scl(1)3381 var % dimids(2) = output_file %dimids_scl(2)3382 var % dimids(3) = output_file %dimids_vel(3)3383 var % dimvarids(1) = output_file %dimvarids_scl(1)3384 var % dimvarids(2) = output_file %dimvarids_scl(2)3385 var % dimvarids(3) = output_file %dimvarids_vel(3)3386 var %to_be_processed = init_variables_required3387 var %is_internal = .FALSE.3388 var %task = "interpolate_3d"3396 var%nt = 1 3397 var%lod = 2 3398 var%ndim = 3 3399 var%dimids(1) = output_file%dimids_scl(1) 3400 var%dimids(2) = output_file%dimids_scl(2) 3401 var%dimids(3) = output_file%dimids_vel(3) 3402 var%dimvarids(1) = output_file%dimvarids_scl(1) 3403 var%dimvarids(2) = output_file%dimvarids_scl(2) 3404 var%dimvarids(3) = output_file%dimvarids_vel(3) 3405 var%to_be_processed = init_variables_required 3406 var%is_internal = .FALSE. 3407 var%task = "interpolate_3d" 3389 3408 3390 3409 CASE( 'init scalar profile', 'init u profile', 'init v profile') 3391 var %nt = 13392 var %lod = 13393 var %ndim = 13394 var % dimids(1) = output_file %dimids_scl(3) !z3395 var % dimvarids(1) = output_file %dimvarids_scl(3) !z3396 var %to_be_processed = init_variables_required3397 var %is_internal = .FALSE.3398 var %task = "average profile"3410 var%nt = 1 3411 var%lod = 1 3412 var%ndim = 1 3413 var%dimids(1) = output_file%dimids_scl(3) !z 3414 var%dimvarids(1) = output_file%dimvarids_scl(3) !z 3415 var%to_be_processed = init_variables_required 3416 var%is_internal = .FALSE. 3417 var%task = "average profile" 3399 3418 3400 3419 CASE( 'init w profile') 3401 var %nt = 13402 var %lod = 13403 var %ndim = 13404 var % dimids(1) = output_file %dimids_vel(3) !z3405 var % dimvarids(1) = output_file %dimvarids_vel(3) !z3406 var %to_be_processed = init_variables_required3407 var %is_internal = .FALSE.3408 var %task = "average profile"3420 var%nt = 1 3421 var%lod = 1 3422 var%ndim = 1 3423 var%dimids(1) = output_file%dimids_vel(3) !z 3424 var%dimvarids(1) = output_file%dimvarids_vel(3) !z 3425 var%to_be_processed = init_variables_required 3426 var%is_internal = .FALSE. 3427 var%task = "average profile" 3409 3428 3410 3429 CASE( 'surface forcing' ) 3411 var %lod = -13412 var %ndim = 33413 var % dimids(3) = output_file %dimid_time3414 var % dimids(1:2) = output_file %dimids_soil(1:2)3415 var % dimvarids(3) = output_file %dimvarid_time3416 var % dimvarids(1:2) = output_file %dimvarids_soil(1:2)3417 var %to_be_processed = surface_forcing_required3418 var %is_internal = .FALSE.3419 var %task = "interpolate_2d"3430 var%lod = -1 3431 var%ndim = 3 3432 var%dimids(3) = output_file%dimid_time 3433 var%dimids(1:2) = output_file%dimids_soil(1:2) 3434 var%dimvarids(3) = output_file%dimvarid_time 3435 var%dimvarids(1:2) = output_file%dimvarids_soil(1:2) 3436 var%to_be_processed = surface_forcing_required 3437 var%is_internal = .FALSE. 3438 var%task = "interpolate_2d" 3420 3439 3421 3440 CASE( 'left scalar', 'right scalar') 3422 var %lod = -13423 var %ndim = 33424 var % dimids(3) = output_file %dimid_time3425 var % dimids(1) = output_file %dimids_scl(2)3426 var % dimids(2) = output_file %dimids_scl(3)3427 var % dimvarids(3) = output_file %dimvarid_time3428 var % dimvarids(1) = output_file %dimvarids_scl(2)3429 var % dimvarids(2) = output_file %dimvarids_scl(3)3430 var %to_be_processed = boundary_variables_required3431 var %is_internal = .FALSE.3432 var %task = "interpolate_3d"3441 var%lod = -1 3442 var%ndim = 3 3443 var%dimids(3) = output_file%dimid_time 3444 var%dimids(1) = output_file%dimids_scl(2) 3445 var%dimids(2) = output_file%dimids_scl(3) 3446 var%dimvarids(3) = output_file%dimvarid_time 3447 var%dimvarids(1) = output_file%dimvarids_scl(2) 3448 var%dimvarids(2) = output_file%dimvarids_scl(3) 3449 var%to_be_processed = boundary_variables_required 3450 var%is_internal = .FALSE. 3451 var%task = "interpolate_3d" 3433 3452 3434 3453 CASE( 'north scalar', 'south scalar') 3435 var %lod = -13436 var %ndim = 33437 var % dimids(3) = output_file %dimid_time3438 var % dimids(1) = output_file %dimids_scl(1)3439 var % dimids(2) = output_file %dimids_scl(3)3440 var % dimvarids(3) = output_file %dimvarid_time3441 var % dimvarids(1) = output_file %dimvarids_scl(1)3442 var % dimvarids(2) = output_file %dimvarids_scl(3)3443 var %to_be_processed = boundary_variables_required3444 var %is_internal = .FALSE.3445 var %task = "interpolate_3d"3454 var%lod = -1 3455 var%ndim = 3 3456 var%dimids(3) = output_file%dimid_time 3457 var%dimids(1) = output_file%dimids_scl(1) 3458 var%dimids(2) = output_file%dimids_scl(3) 3459 var%dimvarids(3) = output_file%dimvarid_time 3460 var%dimvarids(1) = output_file%dimvarids_scl(1) 3461 var%dimvarids(2) = output_file%dimvarids_scl(3) 3462 var%to_be_processed = boundary_variables_required 3463 var%is_internal = .FALSE. 3464 var%task = "interpolate_3d" 3446 3465 3447 3466 CASE( 'top scalar', 'top w' ) 3448 var %lod = -13449 var %ndim = 33450 var % dimids(3) = output_file %dimid_time3451 var % dimids(1) = output_file %dimids_scl(1)3452 var % dimids(2) = output_file %dimids_scl(2)3453 var % dimvarids(3) = output_file %dimvarid_time3454 var % dimvarids(1) = output_file %dimvarids_scl(1)3455 var % dimvarids(2) = output_file %dimvarids_scl(2)3456 var %to_be_processed = boundary_variables_required3457 var %is_internal = .FALSE.3458 var %task = "interpolate_3d"3467 var%lod = -1 3468 var%ndim = 3 3469 var%dimids(3) = output_file%dimid_time 3470 var%dimids(1) = output_file%dimids_scl(1) 3471 var%dimids(2) = output_file%dimids_scl(2) 3472 var%dimvarids(3) = output_file%dimvarid_time 3473 var%dimvarids(1) = output_file%dimvarids_scl(1) 3474 var%dimvarids(2) = output_file%dimvarids_scl(2) 3475 var%to_be_processed = boundary_variables_required 3476 var%is_internal = .FALSE. 3477 var%task = "interpolate_3d" 3459 3478 3460 3479 CASE( 'left u', 'right u' ) 3461 var %lod = -13462 var %ndim = 33463 var % dimids(3) = output_file %dimid_time3464 var % dimids(1) = output_file %dimids_scl(2)3465 var % dimids(2) = output_file %dimids_scl(3)3466 var % dimvarids(3) = output_file %dimvarid_time3467 var % dimvarids(1) = output_file %dimvarids_scl(2)3468 var % dimvarids(2) = output_file %dimvarids_scl(3)3469 var %to_be_processed = boundary_variables_required3470 var %is_internal = .FALSE.3471 var %task = "interpolate_3d"3480 var%lod = -1 3481 var%ndim = 3 3482 var%dimids(3) = output_file%dimid_time 3483 var%dimids(1) = output_file%dimids_scl(2) 3484 var%dimids(2) = output_file%dimids_scl(3) 3485 var%dimvarids(3) = output_file%dimvarid_time 3486 var%dimvarids(1) = output_file%dimvarids_scl(2) 3487 var%dimvarids(2) = output_file%dimvarids_scl(3) 3488 var%to_be_processed = boundary_variables_required 3489 var%is_internal = .FALSE. 3490 var%task = "interpolate_3d" 3472 3491 3473 3492 CASE( 'north u', 'south u' ) 3474 var %lod = -13475 var %ndim = 33476 var % dimids(3) = output_file %dimid_time !t3477 var % dimids(1) = output_file %dimids_vel(1) !x3478 var % dimids(2) = output_file %dimids_scl(3) !z3479 var % dimvarids(3) = output_file %dimvarid_time3480 var % dimvarids(1) = output_file %dimvarids_vel(1)3481 var % dimvarids(2) = output_file %dimvarids_scl(3)3482 var %to_be_processed = boundary_variables_required3483 var %is_internal = .FALSE.3484 var %task = "interpolate_3d"3493 var%lod = -1 3494 var%ndim = 3 3495 var%dimids(3) = output_file%dimid_time !t 3496 var%dimids(1) = output_file%dimids_vel(1) !x 3497 var%dimids(2) = output_file%dimids_scl(3) !z 3498 var%dimvarids(3) = output_file%dimvarid_time 3499 var%dimvarids(1) = output_file%dimvarids_vel(1) 3500 var%dimvarids(2) = output_file%dimvarids_scl(3) 3501 var%to_be_processed = boundary_variables_required 3502 var%is_internal = .FALSE. 3503 var%task = "interpolate_3d" 3485 3504 3486 3505 CASE( 'top u' ) 3487 var %lod = -13488 var %ndim = 33489 var % dimids(3) = output_file %dimid_time !t3490 var % dimids(1) = output_file %dimids_vel(1) !x3491 var % dimids(2) = output_file %dimids_scl(2) !z3492 var % dimvarids(3) = output_file %dimvarid_time3493 var % dimvarids(1) = output_file %dimvarids_vel(1)3494 var % dimvarids(2) = output_file %dimvarids_scl(2)3495 var %to_be_processed = boundary_variables_required3496 var %is_internal = .FALSE.3497 var %task = "interpolate_3d"3506 var%lod = -1 3507 var%ndim = 3 3508 var%dimids(3) = output_file%dimid_time !t 3509 var%dimids(1) = output_file%dimids_vel(1) !x 3510 var%dimids(2) = output_file%dimids_scl(2) !z 3511 var%dimvarids(3) = output_file%dimvarid_time 3512 var%dimvarids(1) = output_file%dimvarids_vel(1) 3513 var%dimvarids(2) = output_file%dimvarids_scl(2) 3514 var%to_be_processed = boundary_variables_required 3515 var%is_internal = .FALSE. 3516 var%task = "interpolate_3d" 3498 3517 3499 3518 CASE( 'left v', 'right v' ) 3500 var %lod = -13501 var %ndim = 33502 var % dimids(3) = output_file %dimid_time3503 var % dimids(1) = output_file %dimids_vel(2)3504 var % dimids(2) = output_file %dimids_scl(3)3505 var % dimvarids(3) = output_file %dimvarid_time3506 var % dimvarids(1) = output_file %dimvarids_vel(2)3507 var % dimvarids(2) = output_file %dimvarids_scl(3)3508 var %to_be_processed = boundary_variables_required3509 var %is_internal = .FALSE.3510 var %task = "interpolate_3d"3519 var%lod = -1 3520 var%ndim = 3 3521 var%dimids(3) = output_file%dimid_time 3522 var%dimids(1) = output_file%dimids_vel(2) 3523 var%dimids(2) = output_file%dimids_scl(3) 3524 var%dimvarids(3) = output_file%dimvarid_time 3525 var%dimvarids(1) = output_file%dimvarids_vel(2) 3526 var%dimvarids(2) = output_file%dimvarids_scl(3) 3527 var%to_be_processed = boundary_variables_required 3528 var%is_internal = .FALSE. 3529 var%task = "interpolate_3d" 3511 3530 3512 3531 CASE( 'north v', 'south v' ) 3513 var %lod = -13514 var %ndim = 33515 var % dimids(3) = output_file %dimid_time !t3516 var % dimids(1) = output_file %dimids_scl(1) !x3517 var % dimids(2) = output_file %dimids_scl(3) !z3518 var % dimvarids(3) = output_file %dimvarid_time3519 var % dimvarids(1) = output_file %dimvarids_scl(1)3520 var % dimvarids(2) = output_file %dimvarids_scl(3)3521 var %to_be_processed = boundary_variables_required3522 var %is_internal = .FALSE.3523 var %task = "interpolate_3d"3532 var%lod = -1 3533 var%ndim = 3 3534 var%dimids(3) = output_file%dimid_time !t 3535 var%dimids(1) = output_file%dimids_scl(1) !x 3536 var%dimids(2) = output_file%dimids_scl(3) !z 3537 var%dimvarids(3) = output_file%dimvarid_time 3538 var%dimvarids(1) = output_file%dimvarids_scl(1) 3539 var%dimvarids(2) = output_file%dimvarids_scl(3) 3540 var%to_be_processed = boundary_variables_required 3541 var%is_internal = .FALSE. 3542 var%task = "interpolate_3d" 3524 3543 3525 3544 CASE( 'top v' ) 3526 var %lod = -13527 var %ndim = 33528 var % dimids(3) = output_file %dimid_time !t3529 var % dimids(1) = output_file %dimids_scl(1) !x3530 var % dimids(2) = output_file %dimids_vel(2) !z3531 var % dimvarids(3) = output_file %dimvarid_time3532 var % dimvarids(1) = output_file %dimvarids_scl(1)3533 var % dimvarids(2) = output_file %dimvarids_vel(2)3534 var %to_be_processed = boundary_variables_required3535 var %is_internal = .FALSE.3536 var %task = "interpolate_3d"3545 var%lod = -1 3546 var%ndim = 3 3547 var%dimids(3) = output_file%dimid_time !t 3548 var%dimids(1) = output_file%dimids_scl(1) !x 3549 var%dimids(2) = output_file%dimids_vel(2) !z 3550 var%dimvarids(3) = output_file%dimvarid_time 3551 var%dimvarids(1) = output_file%dimvarids_scl(1) 3552 var%dimvarids(2) = output_file%dimvarids_vel(2) 3553 var%to_be_processed = boundary_variables_required 3554 var%is_internal = .FALSE. 3555 var%task = "interpolate_3d" 3537 3556 3538 3557 CASE( 'left w', 'right w') 3539 var %lod = -13540 var %ndim = 33541 var % dimids(3) = output_file %dimid_time3542 var % dimids(1) = output_file %dimids_scl(2)3543 var % dimids(2) = output_file %dimids_vel(3)3544 var % dimvarids(3) = output_file %dimvarid_time3545 var % dimvarids(1) = output_file %dimvarids_scl(2)3546 var % dimvarids(2) = output_file %dimvarids_vel(3)3547 var %to_be_processed = boundary_variables_required3548 var %is_internal = .FALSE.3549 var %task = "interpolate_3d"3558 var%lod = -1 3559 var%ndim = 3 3560 var%dimids(3) = output_file%dimid_time 3561 var%dimids(1) = output_file%dimids_scl(2) 3562 var%dimids(2) = output_file%dimids_vel(3) 3563 var%dimvarids(3) = output_file%dimvarid_time 3564 var%dimvarids(1) = output_file%dimvarids_scl(2) 3565 var%dimvarids(2) = output_file%dimvarids_vel(3) 3566 var%to_be_processed = boundary_variables_required 3567 var%is_internal = .FALSE. 3568 var%task = "interpolate_3d" 3550 3569 3551 3570 CASE( 'north w', 'south w' ) 3552 var %lod = -13553 var %ndim = 33554 var % dimids(3) = output_file %dimid_time !t3555 var % dimids(1) = output_file %dimids_scl(1) !x3556 var % dimids(2) = output_file %dimids_vel(3) !z3557 var % dimvarids(3) = output_file %dimvarid_time3558 var % dimvarids(1) = output_file %dimvarids_scl(1)3559 var % dimvarids(2) = output_file %dimvarids_vel(3)3560 var %to_be_processed = boundary_variables_required3561 var %is_internal = .FALSE.3562 var %task = "interpolate_3d"3571 var%lod = -1 3572 var%ndim = 3 3573 var%dimids(3) = output_file%dimid_time !t 3574 var%dimids(1) = output_file%dimids_scl(1) !x 3575 var%dimids(2) = output_file%dimids_vel(3) !z 3576 var%dimvarids(3) = output_file%dimvarid_time 3577 var%dimvarids(1) = output_file%dimvarids_scl(1) 3578 var%dimvarids(2) = output_file%dimvarids_vel(3) 3579 var%to_be_processed = boundary_variables_required 3580 var%is_internal = .FALSE. 3581 var%task = "interpolate_3d" 3563 3582 3564 3583 CASE( 'time series' ) 3565 var %lod = 03566 var %ndim = 13567 var % dimids(1) = output_file %dimid_time !t3568 var % dimvarids(1) = output_file %dimvarid_time3569 var %to_be_processed = .TRUE.3570 var %is_internal = .FALSE.3571 var %task = "average profile"3584 var%lod = 0 3585 var%ndim = 1 3586 var%dimids(1) = output_file%dimid_time !t 3587 var%dimvarids(1) = output_file%dimvarid_time 3588 var%to_be_processed = .TRUE. 3589 var%is_internal = .FALSE. 3590 var%task = "average profile" 3572 3591 3573 3592 CASE( 'constant scalar profile' ) 3574 var %lod = -13575 var %ndim = 23576 var % dimids(2) = output_file %dimid_time !t3577 var % dimids(1) = output_file %dimids_scl(3) !z3578 var % dimvarids(2) = output_file %dimvarid_time3579 var % dimvarids(1) = output_file %dimvarids_scl(3)3580 var %to_be_processed = .TRUE.3581 var %is_internal = .FALSE.3582 var %task = "set profile"3593 var%lod = -1 3594 var%ndim = 2 3595 var%dimids(2) = output_file%dimid_time !t 3596 var%dimids(1) = output_file%dimids_scl(3) !z 3597 var%dimvarids(2) = output_file%dimvarid_time 3598 var%dimvarids(1) = output_file%dimvarids_scl(3) 3599 var%to_be_processed = .TRUE. 3600 var%is_internal = .FALSE. 3601 var%task = "set profile" 3583 3602 3584 3603 CASE( 'large-scale scalar forcing' ) 3585 var %lod = -13586 var %ndim = 23587 var % dimids(2) = output_file %dimid_time !t3588 var % dimids(1) = output_file %dimids_scl(3) !z3589 var % dimvarids(2) = output_file %dimvarid_time3590 var % dimvarids(1) = output_file %dimvarids_scl(3)3591 var %to_be_processed = ls_forcing_variables_required3592 var %is_internal = .FALSE.3593 var %task = "average large-scale profile"3604 var%lod = -1 3605 var%ndim = 2 3606 var%dimids(2) = output_file%dimid_time !t 3607 var%dimids(1) = output_file%dimids_scl(3) !z 3608 var%dimvarids(2) = output_file%dimvarid_time 3609 var%dimvarids(1) = output_file%dimvarids_scl(3) 3610 var%to_be_processed = ls_forcing_variables_required 3611 var%is_internal = .FALSE. 3612 var%task = "average large-scale profile" 3594 3613 3595 3614 CASE( 'geostrophic' ) 3596 var %lod = -13597 var %ndim = 23598 var % dimids(2) = output_file %dimid_time !t3599 var % dimids(1) = output_file %dimids_scl(3) !z3600 var % dimvarids(2) = output_file %dimvarid_time3601 var % dimvarids(1) = output_file %dimvarids_scl(3)3602 var %to_be_processed = .TRUE.3603 var %is_internal = .FALSE.3604 var %task = "geostrophic winds"3615 var%lod = -1 3616 var%ndim = 2 3617 var%dimids(2) = output_file%dimid_time !t 3618 var%dimids(1) = output_file%dimids_scl(3) !z 3619 var%dimvarids(2) = output_file%dimvarid_time 3620 var%dimvarids(1) = output_file%dimvarids_scl(3) 3621 var%to_be_processed = .TRUE. 3622 var%is_internal = .FALSE. 3623 var%task = "geostrophic winds" 3605 3624 3606 3625 CASE( 'large-scale w forcing' ) 3607 var %lod = -13608 var %ndim = 23609 var % dimids(2) = output_file %dimid_time !t3610 var % dimids(1) = output_file %dimids_vel(3) !z3611 var % dimvarids(2) = output_file %dimvarid_time3612 var % dimvarids(1) = output_file %dimvarids_vel(3)3613 var %to_be_processed = ls_forcing_variables_required3614 var %is_internal = .FALSE.3615 var %task = "average large-scale profile"3626 var%lod = -1 3627 var%ndim = 2 3628 var%dimids(2) = output_file%dimid_time !t 3629 var%dimids(1) = output_file%dimids_vel(3) !z 3630 var%dimvarids(2) = output_file%dimvarid_time 3631 var%dimvarids(1) = output_file%dimvarids_vel(3) 3632 var%to_be_processed = ls_forcing_variables_required 3633 var%is_internal = .FALSE. 3634 var%task = "average large-scale profile" 3616 3635 3617 3636 CASE( 'internal profile' ) 3618 var %lod = -13619 var %ndim = 23620 var % dimids(2) = output_file %dimid_time !t3621 var % dimids(1) = output_file %dimids_scl(3) !z3622 var % dimvarids(2) = output_file %dimvarid_time3623 var % dimvarids(1) = output_file %dimvarids_scl(3)3624 var %to_be_processed = .TRUE.3625 var %is_internal = .TRUE.3626 var %task = "internal profile"3637 var%lod = -1 3638 var%ndim = 2 3639 var%dimids(2) = output_file%dimid_time !t 3640 var%dimids(1) = output_file%dimids_scl(3) !z 3641 var%dimvarids(2) = output_file%dimvarid_time 3642 var%dimvarids(1) = output_file%dimvarids_scl(3) 3643 var%to_be_processed = .TRUE. 3644 var%is_internal = .TRUE. 3645 var%task = "internal profile" 3627 3646 3628 3647 CASE DEFAULT … … 3630 3649 CALL inifor_abort ('init_nc_var', message) 3631 3650 3632 3633 3634 3635 3636 3637 3638 3639 CALL report('fini_variables', 'Deallocating variable table', cfg %debug)3640 3641 3642 3643 3644 3645 3646 3647 CALL report('fini_io_groups', 'Deallocating IO groups', cfg %debug)3648 3649 3650 3651 3652 3653 3654 3655 CALL report('fini_file_lists', 'Deallocating file lists', cfg %debug)3656 3657 3658 3651 END SELECT 3652 3653 END FUNCTION init_nc_var 3654 3655 3656 SUBROUTINE fini_variables() 3657 3658 CALL report('fini_variables', 'Deallocating variable table', cfg%debug) 3659 DEALLOCATE( input_var_table ) 3660 3661 END SUBROUTINE fini_variables 3662 3663 3664 SUBROUTINE fini_io_groups() 3665 3666 CALL report('fini_io_groups', 'Deallocating IO groups', cfg%debug) 3667 DEALLOCATE( io_group_list ) 3668 3669 END SUBROUTINE fini_io_groups 3670 3671 3672 SUBROUTINE fini_file_lists() 3673 3674 CALL report('fini_file_lists', 'Deallocating file lists', cfg%debug) 3675 DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files ) 3676 3677 END SUBROUTINE fini_file_lists 3659 3678 3660 3679 … … 3670 3689 !> array will match a COSMO-DE scalar array. 3671 3690 !------------------------------------------------------------------------------! 3672 3673 3674 3675 3676 3677 3678 3679 REAL(dp), ALLOCATABLE :: basic_state_pressure(:)3680 3681 3682 3683 3684 3685 input_buffer(:) %is_preprocessed = .FALSE.3686 3687 SELECT CASE( group %kind )3691 SUBROUTINE preprocess(group, input_buffer, cosmo_grid, iter) 3692 3693 TYPE(io_group), INTENT(INOUT), TARGET :: group 3694 TYPE(container), INTENT(INOUT), ALLOCATABLE :: input_buffer(:) 3695 TYPE(grid_definition), INTENT(IN) :: cosmo_grid 3696 INTEGER, INTENT(IN) :: iter 3697 3698 REAL(wp), ALLOCATABLE :: basic_state_pressure(:) 3699 TYPE(container), ALLOCATABLE :: preprocess_buffer(:) 3700 INTEGER :: hour, dt 3701 INTEGER :: i, j, k 3702 INTEGER :: nx, ny, nz 3703 3704 input_buffer(:)%is_preprocessed = .FALSE. 3705 3706 SELECT CASE( group%kind ) 3688 3707 3689 3708 CASE( 'velocities' ) … … 3694 3713 ! 3695 3714 !-- Allocate u and v arrays with scalar dimensions 3696 nx = SIZE(input_buffer(1) %array, 1)3697 ny = SIZE(input_buffer(1) %array, 2)3698 nz = SIZE(input_buffer(1) %array, 3)3699 ALLOCATE( preprocess_buffer(1) %array(nx, ny, nz) ) ! u buffer3700 ALLOCATE( preprocess_buffer(2) %array(nx, ny, nz) ) ! v buffer3701 3702 CALL run_control('time', 'alloc')3715 nx = SIZE(input_buffer(1)%array, 1) 3716 ny = SIZE(input_buffer(1)%array, 2) 3717 nz = SIZE(input_buffer(1)%array, 3) 3718 ALLOCATE( preprocess_buffer(1)%array(nx, ny, nz) ) ! u buffer 3719 ALLOCATE( preprocess_buffer(2)%array(nx, ny, nz) ) ! v buffer 3720 3721 CALL log_runtime('time', 'alloc') 3703 3722 3704 3723 ! 3705 3724 !-- interpolate U and V to centres 3706 CALL centre_velocities( u_face = input_buffer(1) %array, &3707 v_face = input_buffer(2) %array, &3708 u_centre = preprocess_buffer(1) %array, &3709 v_centre = preprocess_buffer(2) %array )3725 CALL centre_velocities( u_face = input_buffer(1)%array, & 3726 v_face = input_buffer(2)%array, & 3727 u_centre = preprocess_buffer(1)%array, & 3728 v_centre = preprocess_buffer(2)%array ) 3710 3729 3711 cfg %rotation_method = 'rotated-pole'3712 SELECT CASE(cfg %rotation_method)3713 3714 CASE('rotated-pole')3715 ! 3716 !-- rotate U and V to PALM-4U orientation and overwrite U and V with3717 !-- rotated velocities3718 DOk = 1, nz3719 DOj = 1, ny3720 DOi = 1, nx3721 CALL uv2uvrot( urot = preprocess_buffer(1) %array(i,j,k), &3722 vrot = preprocess_buffer(2) %array(i,j,k), &3723 rlat = cosmo_grid %lat(j-1), &3724 rlon = cosmo_grid %lon(i-1), &3725 pollat = phi_cn, &3726 pollon = lambda_cn, &3727 u = input_buffer(1) %array(i,j,k), &3728 v = input_buffer(2) %array(i,j,k) )3729 ENDDO3730 ENDDO3731 ENDDO3732 3733 CASE DEFAULT3734 message = "Rotation method '" // TRIM(cfg %rotation_method) // &3735 "' not recognized."3736 CALL inifor_abort('preprocess', message)3730 cfg%rotation_method = 'rotated-pole' 3731 SELECT CASE(cfg%rotation_method) 3732 3733 CASE('rotated-pole') 3734 ! 3735 !-- rotate U and V to PALM-4U orientation and overwrite U and V with 3736 !-- rotated velocities 3737 DO k = 1, nz 3738 DO j = 1, ny 3739 DO i = 1, nx 3740 CALL uv2uvrot( urot = preprocess_buffer(1)%array(i,j,k), & 3741 vrot = preprocess_buffer(2)%array(i,j,k), & 3742 rlat = cosmo_grid%lat(j-1), & 3743 rlon = cosmo_grid%lon(i-1), & 3744 pollat = phi_cn, & 3745 pollon = lambda_cn, & 3746 u = input_buffer(1)%array(i,j,k), & 3747 v = input_buffer(2)%array(i,j,k) ) 3748 ENDDO 3749 ENDDO 3750 ENDDO 3751 3752 CASE DEFAULT 3753 message = "Rotation method '" // TRIM(cfg%rotation_method) // & 3754 "' not recognized." 3755 CALL inifor_abort('preprocess', message) 3737 3756 3738 3757 END SELECT 3739 3758 3740 input_buffer(1) % array(1,:,:) = 0.0_dp3741 input_buffer(2) % array(1,:,:) = 0.0_dp3742 input_buffer(1) % array(:,1,:) = 0.0_dp3743 input_buffer(2) % array(:,1,:) = 0.0_dp3744 3745 input_buffer(1:2) %is_preprocessed = .TRUE.3746 CALL run_control('time', 'comp')3759 input_buffer(1)%array(1,:,:) = 0.0_wp 3760 input_buffer(2)%array(1,:,:) = 0.0_wp 3761 input_buffer(1)%array(:,1,:) = 0.0_wp 3762 input_buffer(2)%array(:,1,:) = 0.0_wp 3763 3764 input_buffer(1:2)%is_preprocessed = .TRUE. 3765 CALL log_runtime('time', 'comp') 3747 3766 3748 3767 DEALLOCATE( preprocess_buffer ) 3749 CALL run_control('time', 'alloc')3750 3751 message = "Input buffers for group '" // TRIM(group %kind) // "'"//&3768 CALL log_runtime('time', 'alloc') 3769 3770 message = "Input buffers for group '" // TRIM(group%kind) // "'"//& 3752 3771 " preprocessed sucessfully." 3753 3772 CALL report('preprocess', message) 3754 3773 3755 3774 CASE( 'thermodynamics' ) ! T, P, QV 3756 nx = SIZE(input_buffer(1) %array, 1)3757 ny = SIZE(input_buffer(1) %array, 2)3758 nz = SIZE(input_buffer(1) %array, 3)3775 nx = SIZE(input_buffer(1)%array, 1) 3776 ny = SIZE(input_buffer(1)%array, 2) 3777 nz = SIZE(input_buffer(1)%array, 3) 3759 3778 3760 3779 ! 3761 3780 !-- Compute absolute pressure if presure perturbation has been read in. 3762 IF ( TRIM(group % in_var_list(2) %name) == 'PP' ) THEN3781 IF ( TRIM(group%in_var_list(2)%name) == 'PP' ) THEN 3763 3782 message = "Absolute pressure, P, not available, " // & 3764 3783 "computing from pressure preturbation PP." … … 3766 3785 3767 3786 ALLOCATE( basic_state_pressure(1:nz) ) 3768 CALL run_control('time', 'alloc')3769 3770 DO j = 1, ny3771 DO i = 1, nx3772 3773 CALL get_basic_state(cosmo_grid %hfl(i,j,:), BETA, P_SL, T_SL,&3787 CALL log_runtime('time', 'alloc') 3788 3789 DO j = 1, ny 3790 DO i = 1, nx 3791 3792 CALL get_basic_state(cosmo_grid%hfl(i,j,:), BETA, P_SL, T_SL,& 3774 3793 RD, G, basic_state_pressure) 3775 3794 … … 3777 3796 !-- Overwrite pressure perturbation with absolute pressure. HECTO 3778 3797 !-- converts pressure perturbation from hPa to Pa. 3779 input_buffer (2) %array(i,j,:) = &3780 HECTO * input_buffer (2) %array(i,j,:) + &3798 input_buffer (2)%array(i,j,:) = & 3799 HECTO * input_buffer (2)%array(i,j,:) + & 3781 3800 basic_state_pressure(:) 3782 3801 3783 3802 ENDDO 3784 3803 ENDDO 3785 CALL run_control('time', 'comp')3804 CALL log_runtime('time', 'comp') 3786 3805 3787 3806 DEALLOCATE( basic_state_pressure ) 3788 CALL run_control('time', 'alloc')3789 3790 group % in_var_list(2) %name = 'P'3807 CALL log_runtime('time', 'alloc') 3808 3809 group%in_var_list(2)%name = 'P' 3791 3810 3792 3811 ENDIF 3793 3812 ! 3794 3813 !-- mark pressure as preprocessed 3795 input_buffer(2) %is_preprocessed = .TRUE.3814 input_buffer(2)%is_preprocessed = .TRUE. 3796 3815 3797 3816 ! 3798 3817 !-- Copy temperature to the last input buffer array 3799 3818 ALLOCATE( & 3800 input_buffer( group % n_output_quantities ) %array (nx, ny, nz) &3819 input_buffer( group%n_output_quantities )%array (nx, ny, nz) & 3801 3820 ) 3802 input_buffer(group % n_output_quantities) % array(:,:,:) = & 3803 input_buffer(1) % array(:,:,:) 3821 CALL log_runtime('time', 'alloc') 3822 input_buffer(group%n_output_quantities)%array(:,:,:) = & 3823 input_buffer(1)%array(:,:,:) 3804 3824 3805 3825 ! 3806 3826 !-- Convert absolute in place to potential temperature 3807 3827 CALL potential_temperature( & 3808 t = input_buffer(1) %array(:,:,:), &3809 p = input_buffer(2) %array(:,:,:), &3828 t = input_buffer(1)%array(:,:,:), & 3829 p = input_buffer(2)%array(:,:,:), & 3810 3830 p_ref = P_REF, & 3811 3831 r = RD_PALM, & … … 3815 3835 ! 3816 3836 !-- mark potential temperature as preprocessed 3817 input_buffer(1) %is_preprocessed = .TRUE.3837 input_buffer(1)%is_preprocessed = .TRUE. 3818 3838 3819 3839 ! 3820 3840 !-- Convert temperature copy to density 3821 3841 CALL moist_density( & 3822 t_rho = input_buffer(group % n_output_quantities) %array(:,:,:), &3823 p = input_buffer(2) %array(:,:,:), &3824 qv = input_buffer(3) %array(:,:,:), &3842 t_rho = input_buffer(group%n_output_quantities)%array(:,:,:), & 3843 p = input_buffer(2)%array(:,:,:), & 3844 qv = input_buffer(3)%array(:,:,:), & 3825 3845 rd = RD, & 3826 3846 rv = RV & … … 3829 3849 ! 3830 3850 !-- mark qv as preprocessed 3831 input_buffer(3) %is_preprocessed = .TRUE.3851 input_buffer(3)%is_preprocessed = .TRUE. 3832 3852 3833 3853 ! 3834 3854 !-- mark density as preprocessed 3835 input_buffer(group % n_output_quantities) %is_preprocessed = .TRUE.3836 3837 3838 message = "Input buffers for group '" // TRIM(group %kind) // "'"//&3855 input_buffer(group%n_output_quantities)%is_preprocessed = .TRUE. 3856 3857 3858 message = "Input buffers for group '" // TRIM(group%kind) // "'"//& 3839 3859 " preprocessed sucessfully." 3840 3860 CALL report('preprocess', message) 3841 CALL run_control('time', 'comp')3842 3861 3843 3862 CASE( 'scalar' ) ! S or W 3844 input_buffer(:) % is_preprocessed = .TRUE. 3845 CALL run_control('time', 'comp') 3863 input_buffer(:)%is_preprocessed = .TRUE. 3846 3864 3847 3865 CASE( 'soil-temperature' ) ! 3848 3866 3849 CALL fill_water_cells(soiltyp, input_buffer(1) %array, &3850 SIZE(input_buffer(1) %array, 3), &3867 CALL fill_water_cells(soiltyp, input_buffer(1)%array, & 3868 SIZE(input_buffer(1)%array, 3), & 3851 3869 FILL_ITERATIONS) 3852 input_buffer(:) % is_preprocessed = .TRUE. 3853 CALL run_control('time', 'comp') 3870 input_buffer(:)%is_preprocessed = .TRUE. 3854 3871 3855 3872 CASE( 'soil-water' ) ! 3856 3873 3857 CALL fill_water_cells(soiltyp, input_buffer(1) %array, &3858 SIZE(input_buffer(1) %array, 3), &3874 CALL fill_water_cells(soiltyp, input_buffer(1)%array, & 3875 SIZE(input_buffer(1)%array, 3), & 3859 3876 FILL_ITERATIONS) 3860 3877 3861 nx = SIZE(input_buffer(1) %array, 1)3862 ny = SIZE(input_buffer(1) %array, 2)3863 nz = SIZE(input_buffer(1) %array, 3)3864 3865 DO k = 1, nz3866 DO j = 1, ny3867 DO i = 1, nx3868 input_buffer(1) %array(i,j,k) = &3869 input_buffer(1) %array(i,j,k) * d_depth_rho_inv(k)3878 nx = SIZE(input_buffer(1)%array, 1) 3879 ny = SIZE(input_buffer(1)%array, 2) 3880 nz = SIZE(input_buffer(1)%array, 3) 3881 3882 DO k = 1, nz 3883 DO j = 1, ny 3884 DO i = 1, nx 3885 input_buffer(1)%array(i,j,k) = & 3886 input_buffer(1)%array(i,j,k) * d_depth_rho_inv(k) 3870 3887 ENDDO 3871 3888 ENDDO … … 3875 3892 CALL report('preprocess', message) 3876 3893 3877 input_buffer(:) % is_preprocessed = .TRUE. 3878 CALL run_control('time', 'comp') 3894 input_buffer(:)%is_preprocessed = .TRUE. 3879 3895 3880 3896 CASE( 'surface' ) ! 3881 input_buffer(:) % is_preprocessed = .TRUE. 3882 CALL run_control('time', 'comp') 3897 input_buffer(:)%is_preprocessed = .TRUE. 3883 3898 3884 3899 CASE( 'accumulated' ) ! 3885 message = "De-accumulating '" // TRIM(group % in_var_list(1) %name) //&3900 message = "De-accumulating '" // TRIM(group%in_var_list(1)%name) //& 3886 3901 "' in iteration " // TRIM(str(iter)) 3887 3902 CALL report('preprocess', message) … … 3892 3907 3893 3908 ! 3894 !-- input has been accumulated over one hour. Leave as is3895 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral3896 CASE(1)3897 3898 ! 3899 !-- input has been accumulated over two hours. Subtract previous step3900 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral3901 !-- input_buffer(2) %array(:,:,:) carrries two-hour integral3902 CASE(2)3903 CALL deaverage( &3904 avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp, &3905 avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp, &3906 avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )3907 ! 3908 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral of second hour3909 3910 ! 3911 !-- input has been accumulated over three hours. Subtract previous step3912 !-- input_buffer(1) %array(:,:,:) carrries three-hour integral3913 !-- input_buffer(2) %array(:,:,:) still carrries two-hour integral3914 CASE(3)3915 CALL deaverage( &3916 avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp, &3917 avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp, &3918 avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )3919 ! 3920 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral of third hourA3921 3922 CASE DEFAULT3923 message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"3924 CALL inifor_abort('preprocess', message)3909 !-- input has been accumulated over one hour. Leave as is 3910 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral 3911 CASE(1) 3912 3913 ! 3914 !-- input has been accumulated over two hours. Subtract previous step 3915 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral 3916 !-- input_buffer(2)%array(:,:,:) carrries two-hour integral 3917 CASE(2) 3918 CALL deaverage( & 3919 avg_1 = input_buffer(1)%array(:,:,:), t1 = 1.0_wp, & 3920 avg_2 = input_buffer(2)%array(:,:,:), t2 = 1.0_wp, & 3921 avg_3 = input_buffer(1)%array(:,:,:), t3 = 1.0_wp ) 3922 ! 3923 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral of second hour 3924 3925 ! 3926 !-- input has been accumulated over three hours. Subtract previous step 3927 !-- input_buffer(1)%array(:,:,:) carrries three-hour integral 3928 !-- input_buffer(2)%array(:,:,:) still carrries two-hour integral 3929 CASE(3) 3930 CALL deaverage( & 3931 avg_1 = input_buffer(2)%array(:,:,:), t1 = 1.0_wp, & 3932 avg_2 = input_buffer(1)%array(:,:,:), t2 = 1.0_wp, & 3933 avg_3 = input_buffer(1)%array(:,:,:), t3 = 1.0_wp ) 3934 ! 3935 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral of third hourA 3936 3937 CASE DEFAULT 3938 message = "Invalid averaging period '" // TRIM(str(dt)) // " hours" 3939 CALL inifor_abort('preprocess', message) 3925 3940 3926 3941 END SELECT 3927 input_buffer(:) % is_preprocessed = .TRUE. 3928 CALL run_control('time', 'comp') 3942 input_buffer(:)%is_preprocessed = .TRUE. 3929 3943 3930 3944 CASE( 'running average' ) ! 3931 message = "De-averaging '" // TRIM(group % in_var_list(1) %name) // &3945 message = "De-averaging '" // TRIM(group%in_var_list(1)%name) // & 3932 3946 "' in iteration " // TRIM(str(iter)) 3933 3947 CALL report('preprocess', message) … … 3939 3953 SELECT CASE(dt) 3940 3954 ! 3941 !-- input has been accumulated over one hour. Leave as is3942 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral3943 CASE(1)3944 3945 ! 3946 !-- input has been accumulated over two hours. Subtract previous step3947 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral3948 !-- input_buffer(2) %array(:,:,:) carrries two-hour integral3949 CASE(2)3950 CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp, &3951 input_buffer(2) % array(:,:,:), 2.0_dp, &3952 input_buffer(1) % array(:,:,:), 1.0_dp)3953 ! 3954 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral of second hour3955 3956 ! 3957 !-- input has been accumulated over three hours. Subtract previous step3958 !-- input_buffer(1) %array(:,:,:) carrries three-hour integral3959 !-- input_buffer(2) %array(:,:,:) still carrries two-hour integral3960 CASE(3)3961 CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp, &3962 input_buffer(1) % array(:,:,:), 3.0_dp, &3963 input_buffer(1) % array(:,:,:), 1.0_dp)3964 ! 3965 !-- input_buffer(1) %array(:,:,:) carrries one-hour integral of third hourA3966 3967 CASE DEFAULT3968 message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"3969 CALL inifor_abort('preprocess', message)3955 !-- input has been accumulated over one hour. Leave as is 3956 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral 3957 CASE(1) 3958 3959 ! 3960 !-- input has been accumulated over two hours. Subtract previous step 3961 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral 3962 !-- input_buffer(2)%array(:,:,:) carrries two-hour integral 3963 CASE(2) 3964 CALL deaverage( input_buffer(1)%array(:,:,:), 1.0_wp, & 3965 input_buffer(2)%array(:,:,:), 2.0_wp, & 3966 input_buffer(1)%array(:,:,:), 1.0_wp) 3967 ! 3968 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral of second hour 3969 3970 ! 3971 !-- input has been accumulated over three hours. Subtract previous step 3972 !-- input_buffer(1)%array(:,:,:) carrries three-hour integral 3973 !-- input_buffer(2)%array(:,:,:) still carrries two-hour integral 3974 CASE(3) 3975 CALL deaverage( input_buffer(2)%array(:,:,:), 2.0_wp, & 3976 input_buffer(1)%array(:,:,:), 3.0_wp, & 3977 input_buffer(1)%array(:,:,:), 1.0_wp) 3978 ! 3979 !-- input_buffer(1)%array(:,:,:) carrries one-hour integral of third hourA 3980 3981 CASE DEFAULT 3982 message = "Invalid averaging period '" // TRIM(str(dt)) // " hours" 3983 CALL inifor_abort('preprocess', message) 3970 3984 3971 3985 END SELECT 3972 input_buffer(:) %is_preprocessed = .TRUE.3986 input_buffer(:)%is_preprocessed = .TRUE. 3973 3987 3974 3988 CASE DEFAULT 3975 message = "IO group kind '" // TRIM(group %kind) // "' is not supported."3989 message = "IO group kind '" // TRIM(group%kind) // "' is not supported." 3976 3990 CALL inifor_abort('prerpocess', message) 3977 3991 3978 3979 CALL run_control('time', 'comp')3980 3981 3992 END SELECT 3993 CALL log_runtime('time', 'comp') 3994 3995 END SUBROUTINE preprocess 3982 3996 3983 3997 … … 4006 4020 !> array : the soil array (i.e. water content or temperature) 4007 4021 !------------------------------------------------------------------------------! 4008 SUBROUTINE fill_water_cells(soiltyp, array, nz, niter) 4009 INTEGER(hp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp 4010 REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: array 4011 INTEGER, INTENT(IN) :: nz, niter 4012 4013 REAL(dp), DIMENSION(nz) :: column 4014 INTEGER(hp), DIMENSION(:,:), ALLOCATABLE :: old_soiltyp, new_soiltyp 4015 INTEGER :: l, i, j, nx, ny, n_cells, ii, jj, iter 4016 INTEGER, DIMENSION(8) :: di, dj 4017 4018 nx = SIZE(array, 1) 4019 ny = SIZE(array, 2) 4020 di = (/ -1, -1, -1, 0, 0, 1, 1, 1 /) 4021 dj = (/ -1, 0, 1, -1, 1, -1, 0, 1 /) 4022 4023 ALLOCATE(old_soiltyp(SIZE(soiltyp,1), & 4024 SIZE(soiltyp,2) )) 4025 4026 ALLOCATE(new_soiltyp(SIZE(soiltyp,1), & 4027 SIZE(soiltyp,2) )) 4028 4029 old_soiltyp(:,:) = soiltyp(:,:,1) 4030 new_soiltyp(:,:) = soiltyp(:,:,1) 4031 4032 DO iter = 1, niter 4033 4034 DO j = 1, ny 4035 DO i = 1, nx 4036 4037 IF (old_soiltyp(i,j) == WATER_ID) THEN 4038 4039 n_cells = 0 4040 column(:) = 0.0_dp 4041 DO l = 1, SIZE(di) 4042 4043 ii = MIN(nx, MAX(1, i + di(l))) 4044 jj = MIN(ny, MAX(1, j + dj(l))) 4045 4046 IF (old_soiltyp(ii,jj) .NE. WATER_ID) THEN 4047 n_cells = n_cells + 1 4048 column(:) = column(:) + array(ii,jj,:) 4049 ENDIF 4050 4051 ENDDO 4052 4053 ! 4054 !-- Overwrite if at least one non-water neighbour cell is available 4055 IF (n_cells > 0) THEN 4056 array(i,j,:) = column(:) / n_cells 4057 new_soiltyp(i,j) = 0 4022 SUBROUTINE fill_water_cells(soiltyp, array, nz, niter) 4023 INTEGER(iwp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp 4024 REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) :: array 4025 INTEGER, INTENT(IN) :: nz, niter 4026 4027 REAL(wp), DIMENSION(nz) :: column 4028 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: old_soiltyp, new_soiltyp 4029 INTEGER :: l, i, j, nx, ny, n_cells, ii, jj, iter 4030 INTEGER, DIMENSION(8) :: di, dj 4031 4032 nx = SIZE(array, 1) 4033 ny = SIZE(array, 2) 4034 di = (/ -1, -1, -1, 0, 0, 1, 1, 1 /) 4035 dj = (/ -1, 0, 1, -1, 1, -1, 0, 1 /) 4036 4037 ALLOCATE(old_soiltyp(SIZE(soiltyp,1), & 4038 SIZE(soiltyp,2) )) 4039 4040 ALLOCATE(new_soiltyp(SIZE(soiltyp,1), & 4041 SIZE(soiltyp,2) )) 4042 4043 old_soiltyp(:,:) = soiltyp(:,:,1) 4044 new_soiltyp(:,:) = soiltyp(:,:,1) 4045 4046 DO iter = 1, niter 4047 4048 DO j = 1, ny 4049 DO i = 1, nx 4050 4051 IF (old_soiltyp(i,j) == WATER_ID) THEN 4052 4053 n_cells = 0 4054 column(:) = 0.0_wp 4055 DO l = 1, SIZE(di) 4056 4057 ii = MIN(nx, MAX(1, i + di(l))) 4058 jj = MIN(ny, MAX(1, j + dj(l))) 4059 4060 IF (old_soiltyp(ii,jj) .NE. WATER_ID) THEN 4061 n_cells = n_cells + 1 4062 column(:) = column(:) + array(ii,jj,:) 4058 4063 ENDIF 4059 4064 4065 ENDDO 4066 4067 ! 4068 !-- Overwrite if at least one non-water neighbour cell is available 4069 IF (n_cells > 0) THEN 4070 array(i,j,:) = column(:) / n_cells 4071 new_soiltyp(i,j) = 0 4060 4072 ENDIF 4061 4073 4062 ENDDO 4063 ENDDO 4064 4065 old_soiltyp(:,:) = new_soiltyp(:,:) 4074 ENDIF 4066 4075 4067 4076 ENDDO 4068 4069 DEALLOCATE(old_soiltyp, new_soiltyp) 4070 4071 END SUBROUTINE fill_water_cells 4077 ENDDO 4078 4079 old_soiltyp(:,:) = new_soiltyp(:,:) 4080 4081 ENDDO 4082 4083 DEALLOCATE(old_soiltyp, new_soiltyp) 4084 4085 END SUBROUTINE fill_water_cells 4072 4086 4073 4087
Note: See TracChangeset
for help on using the changeset viewer.