!> @file src/grid.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 2017-2018 Leibniz Universitaet Hannover ! Copyright 2017-2018 Deutscher Wetterdienst Offenbach !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: grid.f90 3183 2018-07-27 14:25:55Z raasch $ ! Introduced new PALM grid stretching ! Updated variable names and metadata for PIDS v1.9 compatibility ! Better compatibility with older Intel compilers: ! - avoiding implicit array allocation with new get_netcdf_variable() ! subroutine instead of function ! Improved command line interface: ! - Produce forcing variables when '-mode profile' is being used ! - Renamend initial-condition mode variable 'mode' to 'ic_mode' ! - Improved handling of the start date string ! Removed unnecessary variables and routines ! ! ! 3182 2018-07-27 13:36:03Z suehring ! Initial revision ! ! ! ! Authors: ! -------- ! @author Eckhard Kadasch ! !------------------------------------------------------------------------------! ! Description: ! ------------ !> The grid module contains all variables and routines to specify and work with !> the numerical grids in INIFOR. By convention, all angles are stored in !> radians. !------------------------------------------------------------------------------! MODULE grid USE control USE defs, & ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp, & SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, & BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L USE io, & ONLY: get_netcdf_variable, get_netcdf_attribute, & parse_command_line_arguments, validate_config USE netcdf, & ONLY: NF90_MAX_NAME, NF90_MAX_VAR_DIMS USE transform, & ONLY: rotate_to_cosmo, find_horizontal_neighbours, & compute_horizontal_interp_weights, & find_vertical_neighbours_and_weights, interpolate_2d, & gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d, & project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot USE types USE util IMPLICIT NONE SAVE REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad] REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad] REAL(dp) :: lambda_n = 0.0_dp !< longitude of rotaded pole of COSMO-DE grid [rad] REAL(dp) :: phi_c = 0.0_dp !< rotated-grid latitude of the center of the PALM domain [rad] REAL(dp) :: lambda_c = 0.0_dp !< rotated-grid longitude of the centre of the PALM domain [rad] REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad] REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad] REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m] REAL(dp) :: dy = 0.0_dp !< PALM-4U grid spacing in y direction [m] REAL(dp) :: dz(10) = -1.0_dp !< PALM-4U grid spacing in z direction [m] REAL(dp) :: dz_max = 1000.0_dp !< maximum vertical grid spacing [m] REAL(dp) :: dz_stretch_factor = 1.08_dp !< factor for vertical grid stretching [m] REAL(dp) :: dz_stretch_level = -9999999.9_dp!< height above which the vertical grid will be stretched [m] REAL(dp) :: dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter REAL(dp) :: dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter REAL(dp) :: dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter REAL(dp) :: dxi = 0.0_dp !< inverse PALM-4U grid spacing in x direction [m^-1] REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1] REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1] REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m] REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m] REAL(dp) :: ug = 0.0_dp !< geostrophic wind in x direction [m/s] REAL(dp) :: vg = 0.0_dp !< geostrophic wind in y direction [m/s] REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa] REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m] REAL(dp) :: y0 = 0.0_dp !< y coordinate of PALM-4U Earth tangent [m] REAL(dp) :: z0 = 0.0_dp !< Elevation of the PALM-4U domain above sea level [m] REAL(dp) :: z_top = 0.0_dp !< height of the scalar top boundary [m] REAL(dp) :: zw_top = 0.0_dp !< height of the vertical velocity top boundary [m] REAL(dp) :: lonmin = 0.0_dp !< Minimunm longitude of COSMO-DE's rotated-pole grid REAL(dp) :: lonmax = 0.0_dp !< Maximum longitude of COSMO-DE's rotated-pole grid REAL(dp) :: latmin = 0.0_dp !< Minimunm latitude of COSMO-DE's rotated-pole grid REAL(dp) :: latmax = 0.0_dp !< Maximum latitude of COSMO-DE's rotated-pole grid REAL(dp) :: latitude = 0.0_dp !< geographical latitude of the PALM-4U origin, from inipar namelist [deg] REAL(dp) :: longitude = 0.0_dp !< geographical longitude of the PALM-4U origin, from inipar namelist [deg] REAL(dp) :: origin_lat = 0.0_dp !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg] REAL(dp) :: origin_lon = 0.0_dp !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg] REAL(dp) :: rotation_angle = 0.0_dp !< clockwise angle the PALM-4U north is rotated away from geographical north [deg] REAL(dp) :: end_time = 0.0_dp !< PALM-4U simulation time [s] REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: depths !< COSMO-DE's TERRA-ML soil layer depths REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth_rho_inv !< COSMO-DE's TERRA-ML soil layer thicknesses REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlon !< longitudes of COSMO-DE's rotated-pole grid REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlat !< latitudes of COSMO-DE's rotated-pole grid REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: time !< output times REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: x !< base palm grid x coordinate vector pointed to by grid_definitions REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: xu !< base palm grid xu coordinate vector pointed to by grid_definitions REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: y !< base palm grid y coordinate vector pointed to by grid_definitions REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: yv !< base palm grid yv coordinate vector pointed to by grid_definitions REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z_column !< base palm grid z coordinate vector including the top boundary coordinate (entire column) REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw_column !< base palm grid zw coordinate vector including the top boundary coordinate (entire column) REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z !< base palm grid z coordinate vector pointed to by grid_definitions REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw !< base palm grid zw coordinate vector pointed to by grid_definitions INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: soiltyp !< COSMO-DE soil type map INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched INTEGER :: i !< indexing variable INTEGER :: average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging INTEGER :: k !< indexing variable INTEGER :: nt !< number of output time steps INTEGER :: nx !< number of PALM-4U grid points in x direction INTEGER :: ny !< number of PALM-4U grid points in y direction INTEGER :: nz !< number of PALM-4U grid points in z direction INTEGER :: nlon !< number of longitudal points in target grid (COSMO-DE) INTEGER :: nlat !< number of latitudal points in target grid (COSMO-DE) INTEGER :: nlev !< number of levels in target grid (COSMO-DE) INTEGER :: layers !< number of COSMO-DE soil layers INTEGER :: start_hour_flow !< start of flow forcing in number of hours relative to start_date INTEGER :: start_hour_soil !< start of soil forcing in number of hours relative to start_date, typically equals start_hour_flow INTEGER :: start_hour_radiation !< start of radiation forcing in number of hours relative to start_date, 0 to 2 hours before start_hour_flow to reconstruct hourly averages from one- to three hourly averages of the input data INTEGER :: start_hour_soilmoisture !< start of forcing for the soil moisture spin-up in number of hours relative to start_date, typically -672 (-4 weeks) INTEGER :: end_hour !< simulation time in hours INTEGER :: end_hour_soilmoisture !< end of soil moisture spin-up in hours relative to start_hour_flow INTEGER :: step_hour !< number of hours between forcing time steps LOGICAL :: init_variables_required LOGICAL :: boundary_variables_required LOGICAL :: ls_forcing_variables_required LOGICAL :: profile_grids_required INTEGER :: n_invar = 0 !< number of variables in the input variable table INTEGER :: n_outvar = 0 !< number of variables in the output variable table TYPE(nc_var), ALLOCATABLE, TARGET :: input_var_table(:) !< table of input variables TYPE(nc_var), ALLOCATABLE, TARGET :: output_var_table(:) !< table of input variables TYPE(nc_var) :: cosmo_var !< COSMO-DE dummy variable, used for reading HHL, rlon, rlat TYPE(grid_definition), TARGET :: palm_grid !< PALM-4U grid in the target system (COSMO-DE rotated-pole) TYPE(grid_definition), TARGET :: palm_intermediate !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid TYPE(grid_definition), TARGET :: cosmo_grid !< target system (COSMO-DE rotated-pole) TYPE(grid_definition), TARGET :: scalars_east_grid !< TYPE(grid_definition), TARGET :: scalars_west_grid !< TYPE(grid_definition), TARGET :: scalars_north_grid !< TYPE(grid_definition), TARGET :: scalars_south_grid !< TYPE(grid_definition), TARGET :: scalars_top_grid !< TYPE(grid_definition), TARGET :: scalars_east_intermediate !< TYPE(grid_definition), TARGET :: scalars_west_intermediate !< TYPE(grid_definition), TARGET :: scalars_north_intermediate !< TYPE(grid_definition), TARGET :: scalars_south_intermediate !< TYPE(grid_definition), TARGET :: scalars_top_intermediate !< TYPE(grid_definition), TARGET :: u_initial_grid !< TYPE(grid_definition), TARGET :: u_east_grid !< TYPE(grid_definition), TARGET :: u_west_grid !< TYPE(grid_definition), TARGET :: u_north_grid !< TYPE(grid_definition), TARGET :: u_south_grid !< TYPE(grid_definition), TARGET :: u_top_grid !< TYPE(grid_definition), TARGET :: u_initial_intermediate !< TYPE(grid_definition), TARGET :: u_east_intermediate !< TYPE(grid_definition), TARGET :: u_west_intermediate !< TYPE(grid_definition), TARGET :: u_north_intermediate !< TYPE(grid_definition), TARGET :: u_south_intermediate !< TYPE(grid_definition), TARGET :: u_top_intermediate !< TYPE(grid_definition), TARGET :: v_initial_grid !< TYPE(grid_definition), TARGET :: v_east_grid !< TYPE(grid_definition), TARGET :: v_west_grid !< TYPE(grid_definition), TARGET :: v_north_grid !< TYPE(grid_definition), TARGET :: v_south_grid !< TYPE(grid_definition), TARGET :: v_top_grid !< TYPE(grid_definition), TARGET :: v_initial_intermediate !< TYPE(grid_definition), TARGET :: v_east_intermediate !< TYPE(grid_definition), TARGET :: v_west_intermediate !< TYPE(grid_definition), TARGET :: v_north_intermediate !< TYPE(grid_definition), TARGET :: v_south_intermediate !< TYPE(grid_definition), TARGET :: v_top_intermediate !< TYPE(grid_definition), TARGET :: w_initial_grid !< TYPE(grid_definition), TARGET :: w_east_grid !< TYPE(grid_definition), TARGET :: w_west_grid !< TYPE(grid_definition), TARGET :: w_north_grid !< TYPE(grid_definition), TARGET :: w_south_grid !< TYPE(grid_definition), TARGET :: w_top_grid !< TYPE(grid_definition), TARGET :: w_initial_intermediate !< TYPE(grid_definition), TARGET :: w_east_intermediate !< TYPE(grid_definition), TARGET :: w_west_intermediate !< TYPE(grid_definition), TARGET :: w_north_intermediate !< TYPE(grid_definition), TARGET :: w_south_intermediate !< TYPE(grid_definition), TARGET :: w_top_intermediate !< TYPE(grid_definition), TARGET :: scalar_profile_grid !< TYPE(grid_definition), TARGET :: scalar_profile_intermediate !< TYPE(grid_definition), TARGET :: w_profile_grid !< TYPE(grid_definition), TARGET :: w_profile_intermediate !< TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude, & dz_max, dz_stretch_factor, dz_stretch_level, & !< old grid stretching parameters dz_stretch_level_start, dz_stretch_level_end !< new grid stretching parameters NAMELIST /d3par/ end_time CHARACTER(LEN=LNAME) :: nc_source_text = '' !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...' CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: flow_files CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_moisture_files CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_files CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files CHARACTER(LEN=SNAME) :: input_prefix !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses CHARACTER(LEN=SNAME) :: flow_suffix !< Suffix of flow input files, e.g. 'flow' CHARACTER(LEN=SNAME) :: soil_suffix !< Suffix of soil input files, e.g. 'soil' CHARACTER(LEN=SNAME) :: radiation_suffix !< Suffix of radiation input files, e.g. 'radiation' CHARACTER(LEN=SNAME) :: soilmoisture_suffix !< Suffix of input files for soil moisture spin-up, e.g. 'soilmoisture' TYPE(nc_file) :: output_file TYPE(inifor_config) :: cfg !< Container of INIFOR configuration CONTAINS SUBROUTINE setup_parameters() ! !------------------------------------------------------------------------------ ! Section 1: Define default parameters !------------------------------------------------------------------------------ cfg % start_date = '2013072100' end_hour = 2 start_hour_soil = -2 start_hour_soilmoisture = - (4 * 7 * 24) - 2 lonmin = -5.0_dp * TO_RADIANS lonmax = 5.5_dp * TO_RADIANS latmin = -5.0_dp * TO_RADIANS latmax = 6.5_dp * TO_RADIANS ! COSMO-DE default rotated pole phi_n = 40.0_dp * TO_RADIANS phi_equat = 50.0_dp * TO_RADIANS lambda_n = -170.0_dp * TO_RADIANS ! COMSMO-DE soil layers layers = 8 ALLOCATE( depths(layers), d_depth_rho_inv(layers) ) depths = (/0.005_dp, 0.02_dp, 0.06_dp, 0.18_dp, 0.54_dp, 1.62_dp, 4.86_dp, 14.58_dp/) d_depth_rho_inv = 1.0_dp / & ( (/0.01_dp, 0.02_dp, 0.06_dp, 0.18_dp, 0.54_dp, 1.62_dp, 4.86_dp, 14.58_dp/) * RHO_L ) ! Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g) origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation origin_lon = 13.082744_dp * TO_RADIANS cfg % z0 = 35.0_dp ! Default atmospheric parameters ug = 0.0_dp vg = 0.0_dp cfg % p0 = P_SL ! Parameters for file names start_hour_flow = 0 start_hour_soil = 0 start_hour_radiation = 0 start_hour_soilmoisture = start_hour_flow - 2 end_hour_soilmoisture = start_hour_flow step_hour = FORCING_STEP input_prefix = 'laf' ! 'laf' for COSMO-DE analyses cfg % flow_prefix = input_prefix cfg % soil_prefix = input_prefix cfg % radiation_prefix = input_prefix cfg % soilmoisture_prefix = input_prefix flow_suffix = '-flow' soil_suffix = '-soil' radiation_suffix = '-rad' soilmoisture_suffix = '-soilmoisture' ! !------------------------------------------------------------------------------ ! Section 2: Read command-line arguments, namelist, and grid parameters !------------------------------------------------------------------------------ ! Set default paths and modes cfg % input_path = './' cfg % hhl_file = '' cfg % soiltyp_file = '' cfg % namelist_file = './namelist' cfg % static_driver_file = '' cfg % output_file = './palm-4u-input.nc' cfg % ic_mode = 'volume' cfg % bc_mode = 'real' ! Update default file names and domain centre CALL parse_command_line_arguments( cfg ) output_file % name = cfg % output_file z0 = cfg % z0 p0 = cfg % p0 init_variables_required = .TRUE. boundary_variables_required = TRIM( cfg % bc_mode ) == 'real' ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal' IF ( ls_forcing_variables_required ) THEN message = "Averaging of large-scale forcing profiles " // & "has not been implemented, yet." CALL abort('setup_parameters', message) END IF CALL normalize_path(cfg % input_path) IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc' IF (TRIM(cfg % soiltyp_file) == '') cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc' CALL validate_config( cfg ) CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode)) CALL report('setup_parameters', " forcing mode: " // TRIM(cfg % bc_mode)) CALL report('setup_parameters', " data path: " // TRIM(cfg % input_path)) CALL report('setup_parameters', " hhl file: " // TRIM(cfg % hhl_file)) CALL report('setup_parameters', " soiltyp file: " // TRIM(cfg % soiltyp_file)) CALL report('setup_parameters', " namelist file: " // TRIM(cfg % namelist_file)) CALL report('setup_parameters', " output data file: " // TRIM(output_file % name)) CALL run_control('time', 'init') ! Read in namelist parameters OPEN(10, FILE=cfg % namelist_file) READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz READ(10, NML=d3par) ! end_time CLOSE(10) CALL run_control('time', 'read') end_hour = CEILING( end_time / 3600.0 * step_hour ) ! Generate input file lists CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour, & cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files) CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour, & cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files) CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, & cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files) CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, & cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files) ! !------------------------------------------------------------------------------ ! Section 3: Check for consistency !------------------------------------------------------------------------------ ! !------------------------------------------------------------------------------ ! Section 4: Compute additional parameters !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Section 4.1: COSMO-DE parameters !------------------------------------------------------------------------------ CALL run_control('time', 'init') ! Read COSMO-DE soil type map cosmo_var % name = 'SOILTYP' CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp) message = 'Reading PALM-4U origin from' IF (TRIM(cfg % static_driver_file) .NE. '') THEN origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon') origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat') message = TRIM(message) // " static driver file '" & // TRIM(cfg % static_driver_file) // "'" ELSE origin_lon = longitude origin_lat = latitude message = TRIM(message) // " namlist file '" & // TRIM(cfg % namelist_file) // "'" END IF origin_lon = origin_lon * TO_RADIANS origin_lat = origin_lat * TO_RADIANS CALL report('setup_parameters', message) ! Read in COSMO-DE heights of half layers (vertical cell faces) cosmo_var % name = 'HHL' CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl) CALL run_control('time', 'read') CALL reverse(hhl) nlon = SIZE(hhl, 1) nlat = SIZE(hhl, 2) nlev = SIZE(hhl, 3) CALL run_control('time', 'comp') ! Appoximate COSMO-DE heights of full layers (cell centres) ALLOCATE( hfl(nlon, nlat, nlev-1) ) CALL run_control('time', 'alloc') DO k = 1, nlev-1 hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) + & hhl(:,:,k+1) ) END DO !------------------------------------------------------------------------------ ! Section 4.2: PALM-4U parameters !------------------------------------------------------------------------------ ! PALM-4U domain extents lx = (nx+1) * dx ly = (ny+1) * dy ! PALM-4U point of Earth tangency x0 = 0.0_dp y0 = 0.0_dp ! time vector nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1 ALLOCATE( time(nt) ) CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time) output_file % time => time CALL run_control('time', 'init') ! Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid phi_c = TO_RADIANS * & phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,& phi_n * TO_DEGREES, lambda_n * TO_DEGREES ) lambda_c = TO_RADIANS * & rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,& phi_n * TO_DEGREES, lambda_n * TO_DEGREES, & 0.0_dp ) ! Set gamma according to whether PALM domain is in the northern or southern ! hemisphere of the COSMO-DE rotated-pole system. Gamma assumes either the ! value 0 or PI and is needed to work around around a bug in the rotated-pole ! coordinate transformations. gam = gamma_from_hemisphere(origin_lat, phi_equat) ! Compute the north pole of the rotated-pole grid centred at the PALM-4U domain ! centre. The resulting (phi_cn, lambda_cn) are coordinates in COSMO-DE's ! rotated-pole grid. phi_cn = phic_to_phin(phi_c) lambda_cn = lamc_to_lamn(phi_c, lambda_c) message = "PALM-4U origin:" // NEW_LINE('') // & " lon (lambda) = " // & TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //& " lat (phi ) = " // & TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //& " lon (lambda) = " // & TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & " lat (phi ) = " // & TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" CALL report ('setup_parameters', message) message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // & " lon (lambda) = " // & TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //& " lat (phi ) = " // & TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)" CALL report ('setup_parameters', message) message = "North pole of the rotated palm system:" // NEW_LINE(' ') // & " lon (lambda) = " // & TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // & " lat (phi ) = " // & TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)" CALL report ('setup_parameters', message) CALL run_control('time', 'comp') END SUBROUTINE setup_parameters !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initializes the COSMO-DE, PALM-4U, PALM-4U boundary grids. !------------------------------------------------------------------------------! SUBROUTINE setup_grids() ! setup_grids(inifor_settings(with nx, ny, nz,...)) CHARACTER :: interp_mode !------------------------------------------------------------------------------ ! Section 0: Define base PALM-4U grid coordinate vectors !------------------------------------------------------------------------------ ! palm x y z, we allocate the column to nz+1 in order to include the top ! scalar boundary. The interpolation grids will be associated with ! a shorter column that omits the top element. ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) ) CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x) CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y) CALL stretched_z(z_column, dz, dz_max=dz_max, & dz_stretch_factor=dz_stretch_factor, & dz_stretch_level=dz_stretch_level, & dz_stretch_level_start=dz_stretch_level_start, & dz_stretch_level_end=dz_stretch_level_end, & dz_stretch_factor_array=dz_stretch_factor_array) z(1:nz) = z_column(1:nz) z_top = z_column(nz+1) ! palm xu yv zw, compared to the scalar grid, velocity coordinates ! contain one element less. ALLOCATE( xu(1:nx), yv(1:ny), zw(1:nz-1), zw_column(1:nz)) CALL linspace(dx, lx - dx, xu) CALL linspace(dy, ly - dy, yv) CALL midpoints(z_column, zw_column) zw(1:nz-1) = zw_column(1:nz-1) zw_top = zw_column(nz) !------------------------------------------------------------------------------ ! Section 1: Define initialization and boundary grids !------------------------------------------------------------------------------ CALL init_grid_definition('palm', grid=palm_grid, & xmin=0.0_dp, xmax=lx, & ymin=0.0_dp, ymax=ly, & x0=x0, y0=y0, z0=z0, & nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode) ! Subtracting 1 because arrays will be allocated with nlon + 1 elements. CALL init_grid_definition('cosmo-de', grid=cosmo_grid, & xmin=lonmin, xmax=lonmax, & ymin=latmin, ymax=latmax, & x0=x0, y0=y0, z0=0.0_dp, & nx=nlon-1, ny=nlat-1, nz=nlev-1) ! Define intermediate grid. This is the same as palm_grid except with a much ! coarser vertical grid. The vertical levels are interpolated in each PALM-4U ! column from COSMO-DE's secondary levels. The main levels are then computed as ! the averages of the bounding secondary levels. CALL init_grid_definition('palm intermediate', grid=palm_intermediate, & xmin=0.0_dp, xmax=lx, & ymin=0.0_dp, ymax=ly, & x0=x0, y0=y0, z0=z0, & nx=nx, ny=ny, nz=nlev-2) CALL init_grid_definition('boundary', grid=u_initial_grid, & xmin = dx, xmax = lx - dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = ny, nz = nz, & z=z, ic_mode=cfg % ic_mode) CALL init_grid_definition('boundary', grid=v_initial_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny-1, nz = nz, & z=z, ic_mode=cfg % ic_mode) CALL init_grid_definition('boundary', grid=w_initial_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = nz-1, & z=zw, ic_mode=cfg % ic_mode) CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate, & xmin = dx, xmax = lx - dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny-1, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = nlev - 2) IF (boundary_variables_required) THEN ! !------------------------------------------------------------------------------ ! Section 2: Define PALM-4U boundary grids !------------------------------------------------------------------------------ CALL init_grid_definition('boundary', grid=scalars_east_grid, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz, z=z) CALL init_grid_definition('boundary', grid=scalars_west_grid, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz, z=z) CALL init_grid_definition('boundary', grid=scalars_north_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=scalars_south_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=scalars_top_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = 1, z=(/z_top/)) CALL init_grid_definition('boundary', grid=u_east_grid, & xmin = lx, xmax = lx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz, z=z) CALL init_grid_definition('boundary', grid=u_west_grid, & xmin = 0.0_dp, xmax = 0.0_dp, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz, z=z) CALL init_grid_definition('boundary', grid=u_north_grid, & xmin = dx, xmax = lx - dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=u_south_grid, & xmin = dx, xmax = lx - dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=u_top_grid, & xmin = dx, xmax = lx - dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = ny, nz = 1, z=(/z_top/)) CALL init_grid_definition('boundary', grid=v_east_grid, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny-1, nz = nz, z=z) CALL init_grid_definition('boundary', grid=v_west_grid, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny-1, nz = nz, z=z) CALL init_grid_definition('boundary', grid=v_north_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly, ymax = ly, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=v_south_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.0_dp, ymax = 0.0_dp, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=v_top_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny-1, nz = 1, z=(/z_top/)) CALL init_grid_definition('boundary', grid=w_east_grid, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz - 1, z=zw) CALL init_grid_definition('boundary', grid=w_west_grid, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nz - 1, z=zw) CALL init_grid_definition('boundary', grid=w_north_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz - 1, z=zw) CALL init_grid_definition('boundary', grid=w_south_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nz - 1, z=zw) CALL init_grid_definition('boundary', grid=w_top_grid, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = 1, z=(/zw_top/)) CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate, & xmin = lx, xmax = lx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate, & xmin = 0.0_dp, xmax = 0.0_dp, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate, & xmin = dx, xmax = lx - dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate, & xmin = dx, xmax = lx - dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate, & xmin = dx, xmax = lx - dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx-1, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny-1, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny-1, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly, ymax = ly, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.0_dp, ymax = 0.0_dp, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = dy, ymax = ly - dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny-1, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate, & xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, & xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = ny, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = 0, nz = nlev - 2) CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, & xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, & ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, & x0=x0, y0=y0, z0 = z0, & nx = nx, ny = ny, nz = nlev - 2) END IF ! !------------------------------------------------------------------------------ ! Section 3: Define profile grids !------------------------------------------------------------------------------ profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR. & TRIM(cfg % bc_mode) == 'ideal' ) IF (profile_grids_required) THEN CALL init_grid_definition('boundary', grid=scalar_profile_grid, & xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = 0, nz = nz, z=z) CALL init_grid_definition('boundary', grid=w_profile_grid, & xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = 0, nz = nz - 1, z=zw) CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,& xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = 0, nz = nlev - 2, z=z) CALL init_grid_definition('boundary', grid=w_profile_intermediate, & xmin = 0.5_dp * lx, xmax = 0.5_dp * lx, & ymin = 0.5_dp * ly, ymax = 0.5_dp * ly, & x0=x0, y0=y0, z0 = z0, & nx = 0, ny = 0, nz = nlev - 2, z=zw) END IF ! !------------------------------------------------------------------------------ ! Section 4: Precompute neighbours and weights for interpolation !------------------------------------------------------------------------------ interp_mode = 's' CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode) IF (boundary_variables_required) THEN CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode) END IF interp_mode = 'u' CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) IF (boundary_variables_required) THEN CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode) END IF interp_mode = 'v' CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) IF (boundary_variables_required) THEN CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode) END IF interp_mode = 'w' CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode) IF (boundary_variables_required) THEN CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode) CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode) END IF IF (TRIM(cfg % ic_mode) == 'profile') THEN CALL setup_averaging(cosmo_grid, palm_intermediate, & average_imin, average_imax, average_jmin, average_jmax) END IF END SUBROUTINE setup_grids SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode) TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid TYPE(grid_definition), INTENT(INOUT), TARGET :: palm_grid, palm_intermediate CHARACTER, INTENT(IN) :: kind CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode TYPE(grid_definition), POINTER :: grid REAL(dp), DIMENSION(:), POINTER :: lat, lon REAL(dp), DIMENSION(:,:,:), POINTER :: h LOGICAL :: setup_vertical setup_vertical = .TRUE. IF (PRESENT(ic_mode)) THEN IF (TRIM(ic_mode) == 'profile') setup_vertical = .FALSE. END IF !------------------------------------------------------------------------------ ! Section 1: Horizontal interpolation !------------------------------------------------------------------------------ ! Select horizontal coordinates according to kind of points (s/w, u, v) SELECT CASE(kind) CASE('s') ! scalars lat => cosmo_grid % lat lon => cosmo_grid % lon h => cosmo_grid % hfl CASE('w') ! vertical velocity lat => cosmo_grid % lat lon => cosmo_grid % lon h => cosmo_grid % hhl CASE('u') ! x velocity lat => cosmo_grid % lat lon => cosmo_grid % lonu h => cosmo_grid % hfl CASE('v') ! y velocity lat => cosmo_grid % latv lon => cosmo_grid % lon h => cosmo_grid % hfl CASE DEFAULT message = "Interpolation quantity '" // kind // "' is not supported." CALL abort('setup_interpolation', message) END SELECT grid => palm_intermediate CALL find_horizontal_neighbours(lat, lon, & grid % clat, grid % clon, grid % ii, grid % jj) CALL compute_horizontal_interp_weights(lat, lon, & grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz) !------------------------------------------------------------------------------ ! Section 2: Vertical interpolation !------------------------------------------------------------------------------ IF (setup_vertical) THEN ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) grid % h(:,:,:) = - EARTH_RADIUS ! For w points, use hhl, for scalars use hfl ! compute the full heights for the intermediate grids CALL interpolate_2d(cosmo_grid % hfl, grid % h, grid) CALL find_vertical_neighbours_and_weights(palm_grid, grid) END IF END SUBROUTINE setup_interpolation SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax) TYPE(grid_definition), INTENT(IN) :: cosmo_grid, palm_intermediate INTEGER, INTENT(INOUT) :: imin, imax, jmin, jmax TYPE(grid_definition), POINTER :: grid REAL(dp) :: lonmin_pos,lonmax_pos, latmin_pos, latmax_pos REAL(dp) :: cosmo_dxi, cosmo_dyi cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0)) cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0)) ! find horizontal index ranges for profile averaging lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi imin = FLOOR(lonmin_pos) imax = CEILING(lonmax_pos) jmin = FLOOR(latmin_pos) jmax = CEILING(latmax_pos) ! average heights for intermediate scalar and w profile grids grid => scalar_profile_intermediate ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) grid % h(:,:,:) = - EARTH_RADIUS CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax) CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid) grid => w_profile_intermediate ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) grid % h(:,:,:) = - EARTH_RADIUS CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax) CALL find_vertical_neighbours_and_weights(w_profile_grid, grid) END SUBROUTINE setup_averaging !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initializes grid_definition-type variables. !> !> Input parameters: !> ----------------- !> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids !> as well as grids covering the boundary surfaces. Valid kinds are: !> - 'palm' !> - 'cosmo-de' !> - 'eastwest-scalar' !> !> min, max : Domain minima and maxima in x, y, and z direction. Note !> that these values do not necessarily translate to the outmost coordinates !> of the generated grid but rather refer to the extent of the underlying !> PALM-4U computational domain (i.e. the outer cell faces). The coordinates !> of the generated grid will be inferred from this information taking into !> account the initialization mode ic_mode. For example, the coordinates of a !> boundary grid initialized using mode 'eastwest-scalar' will be located in !> planes one half grid point outwards of xmin and xmax. !> !> z0 : Elevation of the PALM-4U domain above sea level [m] !> !> n : Number of grod points in x, y, and z direction !> !> Output parameters: !> ------------------ !> grid : Grid variable to be initialized. !------------------------------------------------------------------------------! SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, & x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode) CHARACTER(LEN=*), INTENT(IN) :: kind CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode INTEGER, INTENT(IN) :: nx, ny, nz REAL(dp), INTENT(IN) :: xmin, xmax, ymin, ymax REAL(dp), INTENT(IN) :: x0, y0, z0 REAL(dp), INTENT(IN), TARGET, OPTIONAL :: z(:) REAL(dp), INTENT(IN), TARGET, OPTIONAL :: zw(:) TYPE(grid_definition), INTENT(INOUT) :: grid grid % nx = nx grid % ny = ny grid % nz = nz grid % lx = xmax - xmin grid % ly = ymax - ymin grid % x0 = x0 grid % y0 = y0 grid % z0 = z0 SELECT CASE( TRIM (kind) ) CASE('boundary') IF (.NOT.PRESENT(z)) THEN message = "z has not been passed but is required for 'boundary' grids" CALL abort('init_grid_definition', message) END IF ALLOCATE( grid % x(0:nx) ) CALL linspace(xmin, xmax, grid % x) ALLOCATE( grid % y(0:ny) ) CALL linspace(ymin, ymax, grid % y) grid % z => z ! Allocate neighbour indices and weights IF (TRIM(ic_mode) .NE. 'profile') THEN ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) grid % kk(:,:,:,:) = -1 ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) ) grid % w_verti(:,:,:,:) = 0.0_dp END IF CASE('boundary intermediate') ALLOCATE( grid % x(0:nx) ) CALL linspace(xmin, xmax, grid % x) ALLOCATE( grid % y(0:ny) ) CALL linspace(ymin, ymax, grid % y) ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) ) CALL rotate_to_cosmo( & phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude phip = phi_cn, lamp = lambda_cn, & phi = grid % clat, & lam = grid % clon, & gam = gam & ) ! Allocate neighbour indices and weights ALLOCATE( grid % ii(0:nx, 0:ny, 4), & grid % jj(0:nx, 0:ny, 4) ) grid % ii(:,:,:) = -1 grid % jj(:,:,:) = -1 ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) ) grid % w_horiz(:,:,:) = 0.0_dp ! This mode initializes a Cartesian PALM-4U grid and adds the ! corresponding latitudes and longitudes of the rotated pole grid. CASE('palm') IF (.NOT.PRESENT(z)) THEN message = "z has not been passed but is required for 'palm' grids" CALL abort('init_grid_definition', message) END IF IF (.NOT.PRESENT(zw)) THEN message = "zw has not been passed but is required for 'palm' grids" CALL abort('init_grid_definition', message) END IF grid % name(1) = 'x and lon' grid % name(2) = 'y and lat' grid % name(3) = 'z' !TODO: Remove use of global dx, dy, dz variables. Consider !TODO: associating global x,y, and z arrays. ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x) CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y) grid % z => z CALL linspace(xmin + dx, xmax - dx, grid % xu) CALL linspace(ymin + dy, ymax - dy, grid % yv) grid % zw => zw grid % depths => depths ! Allocate neighbour indices and weights IF (TRIM(ic_mode) .NE. 'profile') THEN ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) ) grid % kk(:,:,:,:) = -1 ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) ) grid % w_verti(:,:,:,:) = 0.0_dp END IF CASE('palm intermediate') grid % name(1) = 'x and lon' grid % name(2) = 'y and lat' grid % name(3) = 'interpolated hhl or hfl' !TODO: Remove use of global dx, dy, dz variables. Consider !TODO: associating global x,y, and z arrays. ALLOCATE( grid % x(0:nx), grid % y(0:ny) ) ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) ) CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid % x) CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid % y) CALL linspace(xmin + dx, xmax - dx, grid % xu) CALL linspace(ymin + dy, ymax - dy, grid % yv) grid % depths => depths ! Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) ) ALLOCATE( grid % clonu(1:nx, 0:ny), grid % clatu(1:nx, 0:ny) ) ALLOCATE( grid % clonv(0:nx, 1:ny), grid % clatv(0:nx, 1:ny) ) ! Compute rotated-pole coordinates of... ! ... PALM-4U centres CALL rotate_to_cosmo( & phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude phip = phi_cn, lamp = lambda_cn, & phi = grid % clat, & lam = grid % clon, & gam = gam & ) ! ... PALM-4U u winds CALL rotate_to_cosmo( & phir = project( grid % y, y0, EARTH_RADIUS ), & ! = plate-carree latitude lamr = project( grid % xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude phip = phi_cn, lamp = lambda_cn, & phi = grid % clatu, & lam = grid % clonu, & gam = gam & ) ! ... PALM-4U v winds CALL rotate_to_cosmo( & phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude lamr = project( grid % x, x0, EARTH_RADIUS ), & ! = plate-carree longitude phip = phi_cn, lamp = lambda_cn, & phi = grid % clatv, & lam = grid % clonv, & gam = gam & ) ! Allocate neighbour indices and weights ALLOCATE( grid % ii(0:nx, 0:ny, 4), & grid % jj(0:nx, 0:ny, 4) ) grid % ii(:,:,:) = -1 grid % jj(:,:,:) = -1 ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) ) grid % w_horiz(:,:,:) = 0.0_dp CASE('cosmo-de') grid % name(1) = 'rlon' ! of COMSO-DE cell centres (scalars) grid % name(2) = 'rlat' ! of COMSO-DE cell centres (scalars) grid % name(3) = 'height' ALLOCATE( grid % lon(0:nx), grid % lat(0:ny) ) ALLOCATE( grid % lonu(0:nx), grid % latv(0:ny) ) CALL linspace(xmin, xmax, grid % lon) CALL linspace(ymin, ymax, grid % lat) grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx) grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny) ! Point to heights of half levels (hhl) and compute heights of full ! levels (hfl) as arithmetic averages grid % hhl => hhl grid % hfl => hfl grid % depths => depths CASE DEFAULT message = "Grid kind '" // TRIM(kind) // "' is not recognized." CALL abort('init_grid_definition', message) END SELECT END SUBROUTINE init_grid_definition !------------------------------------------------------------------------------! ! Description: ! ------------ !> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see !> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139 !> !> This routine computes the levels of scalar points. The levels of the velocity !> points are then obtained as the midpoints inbetween using the INIFOR routine !> 'modpoints'. !------------------------------------------------------------------------------! SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, & dz_stretch_level_start, dz_stretch_level_end, & dz_stretch_factor_array) REAL(dp), DIMENSION(:), INTENT(INOUT) :: z, dz, dz_stretch_factor_array REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_level_start, dz_stretch_level_end REAL(dp), INTENT(IN) :: dz_max, dz_stretch_factor, dz_stretch_level INTEGER :: number_stretch_level_start !< number of user-specified start levels for stretching INTEGER :: number_stretch_level_end !< number of user-specified end levels for stretching REAL(dp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end REAL(dp) :: dz_level_end, dz_stretched INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched INTEGER :: dz_stretch_level_index = 0 INTEGER :: k, n, number_dz ! !-- Compute height of u-levels from constant grid length and dz stretch factors IF ( dz(1) == -1.0_dp ) THEN message = 'missing dz' CALL abort( 'stretched_z', message) ELSEIF ( dz(1) <= 0.0_dp ) THEN WRITE( message, * ) 'dz=', dz(1),' <= 0.0' CALL abort( 'stretched_z', message) ENDIF ! !-- Initialize dz_stretch_level_start with the value of dz_stretch_level !-- if it was set by the user IF ( dz_stretch_level /= -9999999.9_dp ) THEN dz_stretch_level_start(1) = dz_stretch_level ENDIF ! !-- Determine number of dz values and stretching levels specified by the !-- user to allow right controlling of the stretching mechanism and to !-- perform error checks. The additional requirement that dz /= dz_max !-- for counting number of user-specified dz values is necessary. Otherwise !-- restarts would abort if the old stretching mechanism with dz_stretch_level !-- is used (Attention: The user is not allowed to specify a dz value equal !-- to the default of dz_max = 999.0). number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max ) number_stretch_level_start = COUNT( dz_stretch_level_start /= & -9999999.9_dp ) number_stretch_level_end = COUNT( dz_stretch_level_end /= & 9999999.9_dp ) ! !-- The number of specified end levels +1 has to be the same than the number !-- of specified dz values IF ( number_dz /= number_stretch_level_end + 1 ) THEN WRITE( message, * ) 'The number of values for dz = ', & number_dz, 'has to be the same than ', & 'the number of values for ', & 'dz_stretch_level_end + 1 = ', & number_stretch_level_end+1 CALL abort( 'stretched_z', message) ENDIF ! !-- The number of specified start levels has to be the same or one less than !-- the number of specified dz values IF ( number_dz /= number_stretch_level_start + 1 .AND. & number_dz /= number_stretch_level_start ) THEN WRITE( message, * ) 'The number of values for dz = ', & number_dz, 'has to be the same or one ', & 'more than& the number of values for ', & 'dz_stretch_level_start = ', & number_stretch_level_start CALL abort( 'stretched_z', message) ENDIF !-- The number of specified start levels has to be the same or one more than !-- the number of specified end levels IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. & number_stretch_level_start /= number_stretch_level_end ) THEN WRITE( message, * ) 'The number of values for ', & 'dz_stretch_level_start = ', & dz_stretch_level_start, 'has to be the ',& 'same or one more than& the number of ', & 'values for dz_stretch_level_end = ', & number_stretch_level_end CALL abort( 'stretched_z', message) ENDIF ! !-- Initialize dz for the free atmosphere with the value of dz_max IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND. & number_stretch_level_start /= 0 ) THEN dz(number_stretch_level_start+1) = dz_max ENDIF ! !-- Initialize the stretching factor if (infinitely) stretching in the free !-- atmosphere is desired (dz_stretch_level_end was not specified for the !-- free atmosphere) IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN dz_stretch_factor_array(number_stretch_level_start) = & dz_stretch_factor ENDIF !-- Allocation of arrays for stretching ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) ) ! !-- The stretching region has to be large enough to allow for a smooth !-- transition between two different grid spacings DO n = 1, number_stretch_level_start min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & 4 * MAX( dz(n),dz(n+1) ) ENDDO IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > & dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN !IF ( ANY( min_dz_stretch_level_end > & ! dz_stretch_level_end ) ) THEN message = 'Each dz_stretch_level_end has to be larger ' // & 'than its corresponding value for ' // & 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& 'to allow for smooth grid stretching' CALL abort('stretched_z', message) ENDIF ! !-- Stretching must not be applied within the prandtl_layer !-- (first two grid points). For the default case dz_stretch_level_start !-- is negative. Therefore the absolut value is checked here. IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',& 'larger than ', dz(1) * 1.5 CALL abort( 'stretched_z', message) ENDIF ! !-- The stretching has to start and end on a grid level. Therefore !-- user-specified values have to ''interpolate'' to the next lowest level IF ( number_stretch_level_start /= 0 ) THEN dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & dz(1)/2.0) / dz(1) ) & * dz(1) + dz(1)/2.0 ENDIF IF ( number_stretch_level_start > 1 ) THEN DO n = 2, number_stretch_level_start dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & dz(n) ) * dz(n) ENDDO ENDIF IF ( number_stretch_level_end /= 0 ) THEN DO n = 1, number_stretch_level_end dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & dz(n+1) ) * dz(n+1) ENDDO ENDIF ! !-- Determine stretching factor if necessary IF ( number_stretch_level_end >= 1 ) THEN CALL calculate_stretching_factor( number_stretch_level_end, dz, & dz_stretch_factor, & dz_stretch_factor_array, & dz_stretch_level_end, & dz_stretch_level_start ) ENDIF z(1) = dz(1) * 0.5_dp ! dz_stretch_level_index = n dz_stretched = dz(1) DO k = 2, n IF ( dz_stretch_level <= z(k-1) .AND. dz_stretched < dz_max ) THEN dz_stretched = dz_stretched * dz_stretch_factor dz_stretched = MIN( dz_stretched, dz_max ) IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1 ENDIF z(k) = z(k-1) + dz_stretched ENDDO !-- Determine u and v height levels considering the possibility of grid !-- stretching in several heights. n = 1 dz_stretch_level_start_index(:) = UBOUND(z, 1) dz_stretch_level_end_index(:) = UBOUND(z, 1) dz_stretched = dz(1) !-- The default value of dz_stretch_level_start is negative, thus the first !-- condition is always true. Hence, the second condition is necessary. DO k = 2, UBOUND(z, 1) IF ( dz_stretch_level_start(n) <= z(k-1) .AND. & dz_stretch_level_start(n) /= -9999999.9_dp ) THEN dz_stretched = dz_stretched * dz_stretch_factor_array(n) IF ( dz(n) > dz(n+1) ) THEN dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz ELSE dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz ENDIF IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) ) & dz_stretch_level_start_index(n) = k-1 ENDIF z(k) = z(k-1) + dz_stretched ! !-- Make sure that the stretching ends exactly at dz_stretch_level_end dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) IF ( dz_level_end < dz(n+1)/3.0 ) THEN z(k) = dz_stretch_level_end(n) dz_stretched = dz(n+1) dz_stretch_level_end_index(n) = k n = n + 1 ENDIF ENDDO DEALLOCATE( min_dz_stretch_level_end ) END SUBROUTINE stretched_z ! Description: [PALM subroutine] ! -----------------------------------------------------------------------------! !> Calculation of the stretching factor through an iterative method. Ideas were !> taken from the paper "Regional stretched grid generation and its application !> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the !> system of equations has two variables (r,l) but four requirements !> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which !> results into an overdetermined system. !------------------------------------------------------------------------------! SUBROUTINE calculate_stretching_factor( number_end, dz, dz_stretch_factor, & dz_stretch_factor_array, & dz_stretch_level_end, & dz_stretch_level_start ) REAL(dp), DIMENSION(:), INTENT(IN) :: dz REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_factor_array REAL(dp), DIMENSION(:), INTENT(IN) :: dz_stretch_level_end, dz_stretch_level_start REAL(dp) :: dz_stretch_factor INTEGER :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached INTEGER :: l_rounded !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2 INTEGER :: n !< loop variable for stretching INTEGER, INTENT(IN) :: number_end !< number of user-specified end levels for stretching REAL(dp) :: delta_l !< absolute difference between l and l_rounded REAL(dp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2 REAL(dp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) REAL(dp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration REAL(dp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region) REAL(dp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly REAL(dp) :: numerator !< numerator of the quotient REAL(dp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly REAL(dp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly REAL(dp) :: dz_stretch_factor_array_2(9) = 1.08_dp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 REAL(dp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors REAL(dp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor REAL(dp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor l = 0 DO n = 1, number_end iterations = 1 stretch_factor_1 = 1.0 stretch_factor_2 = 1.0 delta_total_old = 1.0 IF ( dz(n) > dz(n+1) ) THEN DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) stretch_factor_1 = 1.0 - iterations * stretch_factor_interval distance = ABS( dz_stretch_level_end(n) - & dz_stretch_level_start(n) ) numerator = distance*stretch_factor_1/dz(n) + & stretch_factor_1 - distance/dz(n) IF ( numerator > 0.0 ) THEN l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 l_rounded = NINT( l ) delta_l = ABS( l_rounded - l ) / l ENDIF stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) delta_stretch_factor = ABS( stretch_factor_1 - & stretch_factor_2 ) / & stretch_factor_2 delta_total_new = delta_l + delta_stretch_factor ! !-- stretch_factor_1 is taken to guarantee that the stretching !-- procedure ends as close as possible to dz_stretch_level_end. !-- stretch_factor_2 would guarantee that the stretched dz(n) is !-- equal to dz(n+1) after l_rounded grid levels. IF (delta_total_new < delta_total_old) THEN dz_stretch_factor_array(n) = stretch_factor_1 dz_stretch_factor_array_2(n) = stretch_factor_2 delta_total_old = delta_total_new ENDIF iterations = iterations + 1 ENDDO ELSEIF ( dz(n) < dz(n+1) ) THEN DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) stretch_factor_1 = 1.0 + iterations * stretch_factor_interval distance = ABS( dz_stretch_level_end(n) - & dz_stretch_level_start(n) ) numerator = distance*stretch_factor_1/dz(n) + & stretch_factor_1 - distance/dz(n) l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 l_rounded = NINT( l ) delta_l = ABS( l_rounded - l ) / l stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) delta_stretch_factor = ABS( stretch_factor_1 - & stretch_factor_2 ) / & stretch_factor_2 delta_total_new = delta_l + delta_stretch_factor ! !-- stretch_factor_1 is taken to guarantee that the stretching !-- procedure ends as close as possible to dz_stretch_level_end. !-- stretch_factor_2 would guarantee that the stretched dz(n) is !-- equal to dz(n+1) after l_rounded grid levels. IF (delta_total_new < delta_total_old) THEN dz_stretch_factor_array(n) = stretch_factor_1 dz_stretch_factor_array_2(n) = stretch_factor_2 delta_total_old = delta_total_new ENDIF iterations = iterations + 1 ENDDO ELSE message = 'Two adjacent values of dz must be different' CALL abort( 'calculate_stretching_factor', message) ENDIF ! !-- Check if also the second stretching factor fits into the allowed !-- interval. If not, print a warning for the user. IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR. & dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN WRITE( message, * ) 'stretch_factor_2 = ', & dz_stretch_factor_array_2(n), ' which is',& ' responsible for exactly reaching& dz =',& dz(n+1), 'after a specific amount of', & ' grid levels& exceeds the upper', & ' limit =', stretch_factor_upper_limit, & ' &or lower limit = ', & stretch_factor_lower_limit CALL abort( 'calculate_stretching_factor', message ) ENDIF ENDDO END SUBROUTINE calculate_stretching_factor SUBROUTINE midpoints(z, zw) REAL(dp), INTENT(IN) :: z(0:) REAL(dp), INTENT(OUT) :: zw(1:) INTEGER :: k DO k = 1, UBOUND(zw, 1) zw(k) = 0.5_dp * (z(k-1) + z(k)) END DO END SUBROUTINE midpoints SUBROUTINE setup_io_groups() INTEGER :: ngroups ngroups = 16 ALLOCATE( io_group_list(ngroups) ) !soil temp io_group_list(1) = init_io_group( & in_files = soil_files, & out_vars = output_var_table(1:1), & in_var_list = input_var_table(1:1), & kind = 'soil-temperature' & ) !soil water io_group_list(2) = init_io_group( & in_files = soil_files, & out_vars = output_var_table(2:2), & in_var_list = input_var_table(2:2), & kind = 'soil-water' & ) !potential temperature, surface pressure, including nudging and subsidence io_group_list(3) = init_io_group( & in_files = flow_files, & out_vars = [output_var_table(3:8), output_var_table(42:42), & output_var_table(49:51)], & in_var_list = (/input_var_table(3), input_var_table(17)/), & kind = 'temperature' & ) !specific humidity including nudging and subsidence io_group_list(4) = init_io_group( & in_files = flow_files, & out_vars = [output_var_table(9:14), output_var_table(52:54)], & in_var_list = input_var_table(4:4), & kind = 'scalar' & ) !u and v velocity and geostrophic winds ug, vg io_group_list(5) = init_io_group( & in_files = flow_files, & out_vars = [output_var_table(15:26), output_var_table(43:46)], & !out_vars = output_var_table(15:20), & in_var_list = input_var_table(5:6), & !in_var_list = input_var_table(5:5), & kind = 'velocities' & ) !!v velocity, deprecated! !io_group_list(6) = init_io_group( & ! in_files = flow_files, & ! out_vars = output_var_table(21:26), & ! in_var_list = input_var_table(6:6), & ! kind = 'horizontal velocity' & !) !io_group_list(6) % to_be_processed = .FALSE. !w velocity and subsidence and w nudging io_group_list(7) = init_io_group( & in_files = flow_files, & out_vars = [output_var_table(27:32), output_var_table(47:48)], & in_var_list = input_var_table(7:7), & kind = 'scalar' & ) !rain io_group_list(8) = init_io_group( & in_files = soil_moisture_files, & out_vars = output_var_table(33:33), & in_var_list = input_var_table(8:8), & kind = 'accumulated' & ) io_group_list(8) % to_be_processed = .FALSE. !snow io_group_list(9) = init_io_group( & in_files = soil_moisture_files, & out_vars = output_var_table(34:34), & in_var_list = input_var_table(9:9), & kind = 'accumulated' & ) io_group_list(9) % to_be_processed = .FALSE. !graupel io_group_list(10) = init_io_group( & in_files = soil_moisture_files, & out_vars = output_var_table(35:35), & in_var_list = input_var_table(10:10), & kind = 'accumulated' & ) io_group_list(10) % to_be_processed = .FALSE. !evapotranspiration io_group_list(11) = init_io_group( & in_files = soil_moisture_files, & out_vars = output_var_table(37:37), & in_var_list = input_var_table(11:11), & kind = 'accumulated' & ) io_group_list(11) % to_be_processed = .FALSE. !2m air temperature io_group_list(12) = init_io_group( & in_files = soil_moisture_files, & out_vars = output_var_table(36:36), & in_var_list = input_var_table(12:12), & kind = 'surface' & ) !incoming diffusive sw flux io_group_list(13) = init_io_group( & in_files = radiation_files, & out_vars = output_var_table(38:38), & in_var_list = input_var_table(13:13), & kind = 'running average' & ) io_group_list(13) % to_be_processed = .FALSE. !incoming direct sw flux io_group_list(14) = init_io_group( & in_files = radiation_files, & out_vars = output_var_table(39:39), & in_var_list = input_var_table(14:14), & kind = 'running average' & ) io_group_list(14) % to_be_processed = .FALSE. !sw radiation balance io_group_list(15) = init_io_group( & in_files = radiation_files, & out_vars = output_var_table(40:40), & in_var_list = input_var_table(15:15), & kind = 'running average' & ) io_group_list(15) % to_be_processed = .FALSE. !lw radiation balance io_group_list(16) = init_io_group( & in_files = radiation_files, & out_vars = output_var_table(41:41), & in_var_list = input_var_table(16:16), & kind = 'running average' & ) io_group_list(16) % to_be_processed = .FALSE. END SUBROUTINE setup_io_groups FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group) CHARACTER(LEN=PATH), INTENT(IN) :: in_files(:) CHARACTER(LEN=*), INTENT(IN) :: kind TYPE(nc_var), INTENT(IN) :: out_vars(:) TYPE(nc_var), INTENT(IN) :: in_var_list(:) TYPE(io_group) :: group group % nt = SIZE(in_files) group % nv = SIZE(out_vars) group % kind = TRIM(kind) ALLOCATE(group % in_var_list( SIZE(in_var_list) )) ALLOCATE(group % in_files(group % nt)) ALLOCATE(group % out_vars(group % nv)) group % in_var_list = in_var_list group % in_files = in_files group % out_vars = out_vars group % to_be_processed = .TRUE. END FUNCTION init_io_group !------------------------------------------------------------------------------! ! Description: ! ------------ !> Deallocates all allocated variables. !------------------------------------------------------------------------------! SUBROUTINE fini_grids() CALL report('fini_grids', 'Deallocating grids') DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column) DEALLOCATE(palm_grid%x, palm_grid%y, palm_grid%z, & palm_grid%xu, palm_grid%yv, palm_grid%zw, & palm_grid%clon, palm_grid%clat, & palm_grid%clonu, palm_grid%clatu) DEALLOCATE(palm_intermediate%x, palm_intermediate%y, palm_intermediate%z, & palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,& palm_intermediate%clon, palm_intermediate%clat, & palm_intermediate%clonu, palm_intermediate%clatu) DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat, & cosmo_grid%lonu, cosmo_grid%latv, & cosmo_grid%hfl) END SUBROUTINE fini_grids !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initializes the the variable list. !------------------------------------------------------------------------------! SUBROUTINE setup_variable_tables(ic_mode) CHARACTER(LEN=*), INTENT(IN) :: ic_mode TYPE(nc_var), POINTER :: var IF (TRIM(cfg % start_date) == '') THEN message = 'Simulation start date has not been set.' CALL abort('setup_variable_tables', message) END IF nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date) n_invar = 17 n_outvar = 55 ALLOCATE( input_var_table(n_invar) ) ALLOCATE( output_var_table(n_outvar) ) ! !------------------------------------------------------------------------------ !- Section 1: NetCDF input variables !------------------------------------------------------------------------------ var => input_var_table(1) var % name = 'T_SO' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(2) var % name = 'W_SO' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(3) var % name = 'T' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. var => input_var_table(4) var % name = 'QV' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. var => input_var_table(5) var % name = 'U' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. var => input_var_table(6) var % name = 'V' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. var => input_var_table(7) var % name = 'W' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. var => input_var_table(8) var % name = 'RAIN_GSP' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(9) var % name = 'SNOW_GSP' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(10) var % name = 'GRAU_GSP' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(11) var % name = 'AEVAP_S' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(12) var % name = 'T_2M' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(13) var % name = 'ASWDIFD_S' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(14) var % name = 'ASWDIR_S' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(15) var % name = 'ASOB_S' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(16) var % name = 'ATHB_S' var % to_be_processed = .TRUE. var % is_upside_down = .FALSE. var => input_var_table(17) var % name = 'P' var % to_be_processed = .TRUE. var % is_upside_down = .TRUE. ! !------------------------------------------------------------------------------ !- Section 2: NetCDF output variables !------------------------------------------------------------------------------ ! !------------------------------------------------------------------------------ ! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions !------------------------------------------------------------------------------ output_var_table(1) = init_nc_var( & name = 'init_soil_t', & std_name = "", & long_name = "initial soil temperature", & units = "K", & kind = "init soil", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(2) = init_nc_var( & name = 'init_soil_m', & std_name = "", & long_name = "initial soil moisture", & units = "m^3/m^3", & kind = "init soil", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(3) = init_nc_var( & name = 'init_atmosphere_pt', & std_name = "", & long_name = "initial potential temperature", & units = "K", & kind = "init scalar", & input_id = 1, & ! first in (T, p) IO group output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate, & is_profile = (TRIM(ic_mode) == 'profile') & ) IF (TRIM(ic_mode) == 'profile') THEN output_var_table(3) % grid => scalar_profile_grid output_var_table(3) % intermediate_grid => scalar_profile_intermediate END IF output_var_table(4) = init_nc_var( & name = 'ls_forcing_left_pt', & std_name = "", & long_name = "large-scale forcing for left model boundary for the potential temperature", & units = "K", & kind = "left scalar", & input_id = 1, & grid = scalars_west_grid, & intermediate_grid = scalars_west_intermediate, & output_file = output_file & ) output_var_table(5) = init_nc_var( & name = 'ls_forcing_right_pt', & std_name = "", & long_name = "large-scale forcing for right model boundary for the potential temperature", & units = "K", & kind = "right scalar", & input_id = 1, & grid = scalars_east_grid, & intermediate_grid = scalars_east_intermediate, & output_file = output_file & ) output_var_table(6) = init_nc_var( & name = 'ls_forcing_north_pt', & std_name = "", & long_name = "large-scale forcing for north model boundary for the potential temperature", & units = "K", & kind = "north scalar", & input_id = 1, & grid = scalars_north_grid, & intermediate_grid = scalars_north_intermediate, & output_file = output_file & ) output_var_table(7) = init_nc_var( & name = 'ls_forcing_south_pt', & std_name = "", & long_name = "large-scale forcing for south model boundary for the potential temperature", & units = "K", & kind = "south scalar", & input_id = 1, & grid = scalars_south_grid, & intermediate_grid = scalars_south_intermediate, & output_file = output_file & ) output_var_table(8) = init_nc_var( & name = 'ls_forcing_top_pt', & std_name = "", & long_name = "large-scale forcing for top model boundary for the potential temperature", & units = "K", & kind = "top scalar", & input_id = 1, & grid = scalars_top_grid, & intermediate_grid = scalars_top_intermediate, & output_file = output_file & ) output_var_table(9) = init_nc_var( & name = 'init_atmosphere_qv', & std_name = "", & long_name = "initial specific humidity", & units = "kg/kg", & kind = "init scalar", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate, & is_profile = (TRIM(ic_mode) == 'profile') & ) IF (TRIM(ic_mode) == 'profile') THEN output_var_table(9) % grid => scalar_profile_grid output_var_table(9) % intermediate_grid => scalar_profile_intermediate END IF output_var_table(10) = init_nc_var( & name = 'ls_forcing_left_qv', & std_name = "", & long_name = "large-scale forcing for left model boundary for the specific humidity", & units = "kg/kg", & kind = "left scalar", & input_id = 1, & output_file = output_file, & grid = scalars_west_grid, & intermediate_grid = scalars_west_intermediate & ) output_var_table(11) = init_nc_var( & name = 'ls_forcing_right_qv', & std_name = "", & long_name = "large-scale forcing for right model boundary for the specific humidity", & units = "kg/kg", & kind = "right scalar", & input_id = 1, & output_file = output_file, & grid = scalars_east_grid, & intermediate_grid = scalars_east_intermediate & ) output_var_table(12) = init_nc_var( & name = 'ls_forcing_north_qv', & std_name = "", & long_name = "large-scale forcing for north model boundary for the specific humidity", & units = "kg/kg", & kind = "north scalar", & input_id = 1, & output_file = output_file, & grid = scalars_north_grid, & intermediate_grid = scalars_north_intermediate & ) output_var_table(13) = init_nc_var( & name = 'ls_forcing_south_qv', & std_name = "", & long_name = "large-scale forcing for south model boundary for the specific humidity", & units = "kg/kg", & kind = "south scalar", & input_id = 1, & output_file = output_file, & grid = scalars_south_grid, & intermediate_grid = scalars_south_intermediate & ) output_var_table(14) = init_nc_var( & name = 'ls_forcing_top_qv', & std_name = "", & long_name = "large-scale forcing for top model boundary for the specific humidity", & units = "kg/kg", & kind = "top scalar", & input_id = 1, & output_file = output_file, & grid = scalars_top_grid, & intermediate_grid = scalars_top_intermediate & ) output_var_table(15) = init_nc_var( & name = 'init_atmosphere_u', & std_name = "", & long_name = "initial wind component in x direction", & units = "m/s", & kind = "init u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_initial_grid, & intermediate_grid = u_initial_intermediate, & is_profile = (TRIM(ic_mode) == 'profile') & ) IF (TRIM(ic_mode) == 'profile') THEN output_var_table(15) % grid => scalar_profile_grid output_var_table(15) % intermediate_grid => scalar_profile_intermediate END IF output_var_table(16) = init_nc_var( & name = 'ls_forcing_left_u', & std_name = "", & long_name = "large-scale forcing for left model boundary for the wind component in x direction", & units = "m/s", & kind = "left u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_west_grid, & intermediate_grid = u_west_intermediate & ) output_var_table(17) = init_nc_var( & name = 'ls_forcing_right_u', & std_name = "", & long_name = "large-scale forcing for right model boundary for the wind component in x direction", & units = "m/s", & kind = "right u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_east_grid, & intermediate_grid = u_east_intermediate & ) output_var_table(18) = init_nc_var( & name = 'ls_forcing_north_u', & std_name = "", & long_name = "large-scale forcing for north model boundary for the wind component in x direction", & units = "m/s", & kind = "north u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_north_grid, & intermediate_grid = u_north_intermediate & ) output_var_table(19) = init_nc_var( & name = 'ls_forcing_south_u', & std_name = "", & long_name = "large-scale forcing for south model boundary for the wind component in x direction", & units = "m/s", & kind = "south u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_south_grid, & intermediate_grid = u_south_intermediate & ) output_var_table(20) = init_nc_var( & name = 'ls_forcing_top_u', & std_name = "", & long_name = "large-scale forcing for top model boundary for the wind component in x direction", & units = "m/s", & kind = "top u", & input_id = 1, & ! first in (U, V) I/O group output_file = output_file, & grid = u_top_grid, & intermediate_grid = u_top_intermediate & ) output_var_table(21) = init_nc_var( & name = 'init_atmosphere_v', & std_name = "", & long_name = "initial wind component in y direction", & units = "m/s", & kind = "init v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_initial_grid, & intermediate_grid = v_initial_intermediate, & is_profile = (TRIM(ic_mode) == 'profile') & ) IF (TRIM(ic_mode) == 'profile') THEN output_var_table(21) % grid => scalar_profile_grid output_var_table(21) % intermediate_grid => scalar_profile_intermediate END IF output_var_table(22) = init_nc_var( & name = 'ls_forcing_left_v', & std_name = "", & long_name = "large-scale forcing for left model boundary for the wind component in y direction", & units = "m/s", & kind = "right v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_west_grid, & intermediate_grid = v_west_intermediate & ) output_var_table(23) = init_nc_var( & name = 'ls_forcing_right_v', & std_name = "", & long_name = "large-scale forcing for right model boundary for the wind component in y direction", & units = "m/s", & kind = "right v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_east_grid, & intermediate_grid = v_east_intermediate & ) output_var_table(24) = init_nc_var( & name = 'ls_forcing_north_v', & std_name = "", & long_name = "large-scale forcing for north model boundary for the wind component in y direction", & units = "m/s", & kind = "north v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_north_grid, & intermediate_grid = v_north_intermediate & ) output_var_table(25) = init_nc_var( & name = 'ls_forcing_south_v', & std_name = "", & long_name = "large-scale forcing for south model boundary for the wind component in y direction", & units = "m/s", & kind = "south v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_south_grid, & intermediate_grid = v_south_intermediate & ) output_var_table(26) = init_nc_var( & name = 'ls_forcing_top_v', & std_name = "", & long_name = "large-scale forcing for top model boundary for the wind component in y direction", & units = "m/s", & kind = "top v", & input_id = 2, & ! second in (U, V) I/O group output_file = output_file, & grid = v_top_grid, & intermediate_grid = v_top_intermediate & ) output_var_table(27) = init_nc_var( & name = 'init_atmosphere_w', & std_name = "", & long_name = "initial wind component in z direction", & units = "m/s", & kind = "init w", & input_id = 1, & output_file = output_file, & grid = w_initial_grid, & intermediate_grid = w_initial_intermediate, & is_profile = (TRIM(ic_mode) == 'profile') & ) IF (TRIM(ic_mode) == 'profile') THEN output_var_table(27) % grid => w_profile_grid output_var_table(27) % intermediate_grid => w_profile_intermediate END IF output_var_table(28) = init_nc_var( & name = 'ls_forcing_left_w', & std_name = "", & long_name = "large-scale forcing for left model boundary for the wind component in z direction", & units = "m/s", & kind = "left w", & input_id = 1, & output_file = output_file, & grid = w_west_grid, & intermediate_grid = w_west_intermediate & ) output_var_table(29) = init_nc_var( & name = 'ls_forcing_right_w', & std_name = "", & long_name = "large-scale forcing for right model boundary for the wind component in z direction", & units = "m/s", & kind = "right w", & input_id = 1, & output_file = output_file, & grid = w_east_grid, & intermediate_grid = w_east_intermediate & ) output_var_table(30) = init_nc_var( & name = 'ls_forcing_north_w', & std_name = "", & long_name = "large-scale forcing for north model boundary for the wind component in z direction", & units = "m/s", & kind = "north w", & input_id = 1, & output_file = output_file, & grid = w_north_grid, & intermediate_grid = w_north_intermediate & ) output_var_table(31) = init_nc_var( & name = 'ls_forcing_south_w', & std_name = "", & long_name = "large-scale forcing for south model boundary for the wind component in z direction", & units = "m/s", & kind = "south w", & input_id = 1, & output_file = output_file, & grid = w_south_grid, & intermediate_grid = w_south_intermediate & ) output_var_table(32) = init_nc_var( & name = 'ls_forcing_top_w', & std_name = "", & long_name = "large-scale forcing for top model boundary for the wind component in z direction", & units = "m/s", & kind = "top w", & input_id = 1, & output_file = output_file, & grid = w_top_grid, & intermediate_grid = w_top_intermediate & ) output_var_table(33) = init_nc_var( & name = 'ls_forcing_soil_rain', & std_name = "", & long_name = "large-scale forcing rain", & units = "kg/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(34) = init_nc_var( & name = 'ls_forcing_soil_snow', & std_name = "", & long_name = "large-scale forcing snow", & units = "kg/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(35) = init_nc_var( & name = 'ls_forcing_soil_graupel', & std_name = "", & long_name = "large-scale forcing graupel", & units = "kg/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(36) = init_nc_var( & name = 'ls_forcing_soil_t_2m', & std_name = "", & long_name = "large-scale forcing 2m air temperature", & units = "kg/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(37) = init_nc_var( & name = 'ls_forcing_soil_evap', & std_name = "", & long_name = "large-scale forcing evapo-transpiration", & units = "kg/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) ! Radiation fluxes and balances output_var_table(38) = init_nc_var( & name = 'rad_swd_dif_0', & std_name = "", & long_name = "incoming diffuse shortwave radiative flux at the surface", & units = "W/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(39) = init_nc_var( & name = 'rad_swd_dir_0', & std_name = "", & long_name = "incoming direct shortwave radiative flux at the surface", & units = "W/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(40) = init_nc_var( & name = 'rad_sw_bal_0', & std_name = "", & long_name = "shortwave radiation balance at the surface", & units = "W/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(41) = init_nc_var( & name = 'rad_lw_bal_0', & std_name = "", & long_name = "longwave radiation balance at the surface", & units = "W/m2", & kind = "surface forcing", & input_id = 1, & output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) ! !------------------------------------------------------------------------------ ! Section 2.2: Idealized large-scale forcings !------------------------------------------------------------------------------ output_var_table(42) = init_nc_var( & name = 'surface_forcing_surface_pressure', & std_name = "", & long_name = "surface pressure", & units = "Pa", & kind = "time series", & input_id = 2, & ! second in (T, p) I/O group output_file = output_file, & grid = palm_grid, & intermediate_grid = palm_intermediate & ) output_var_table(43) = init_nc_var( & name = 'ls_forcing_ug', & std_name = "", & long_name = "geostrophic wind (u component)", & units = "m/s", & kind = "constant scalar profile", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(44) = init_nc_var( & name = 'ls_forcing_vg', & std_name = "", & long_name = "geostrophic wind (v component)", & units = "m/s", & kind = "constant scalar profile", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(45) = init_nc_var( & name = 'nudging_u', & std_name = "", & long_name = "wind component in x direction", & units = "m/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(45) % to_be_processed = ls_forcing_variables_required output_var_table(46) = init_nc_var( & name = 'nudging_v', & std_name = "", & long_name = "wind component in y direction", & units = "m/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(46) % to_be_processed = ls_forcing_variables_required output_var_table(47) = init_nc_var( & name = 'ls_forcing_sub_w', & std_name = "", & long_name = "subsidence velocity of w", & units = "m/s", & kind = "large-scale w forcing", & input_id = 1, & output_file = output_file, & grid = w_profile_grid, & intermediate_grid = w_profile_intermediate & ) output_var_table(47) % to_be_processed = ls_forcing_variables_required output_var_table(48) = init_nc_var( & name = 'nudging_w', & std_name = "", & long_name = "wind component in w direction", & units = "m/s", & kind = "large-scale w forcing", & input_id = 1, & output_file = output_file, & grid = w_profile_grid, & intermediate_grid = w_profile_intermediate & ) output_var_table(48) % to_be_processed = ls_forcing_variables_required output_var_table(49) = init_nc_var( & name = 'ls_forcing_adv_pt', & std_name = "", & long_name = "advection of potential temperature", & units = "K/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(49) % to_be_processed = ls_forcing_variables_required output_var_table(50) = init_nc_var( & name = 'ls_forcing_sub_pt', & std_name = "", & long_name = "subsidence velocity of potential temperature", & units = "K/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(50) % to_be_processed = ls_forcing_variables_required output_var_table(51) = init_nc_var( & name = 'nudging_pt', & std_name = "", & long_name = "potential temperature", & units = "K", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(51) % to_be_processed = ls_forcing_variables_required output_var_table(52) = init_nc_var( & name = 'ls_forcing_adv_qv', & std_name = "", & long_name = "advection of specific humidity", & units = "kg/kg/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(52) % to_be_processed = ls_forcing_variables_required output_var_table(53) = init_nc_var( & name = 'ls_forcing_sub_qv', & std_name = "", & long_name = "subsidence velocity of specific humidity", & units = "kg/kg/s", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(53) % to_be_processed = ls_forcing_variables_required output_var_table(54) = init_nc_var( & name = 'nudging_qv', & std_name = "", & long_name = "specific humidity", & units = "kg/kg", & kind = "large-scale scalar forcing", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(54) % to_be_processed = ls_forcing_variables_required output_var_table(55) = init_nc_var( & name = 'nudging_tau', & std_name = "", & long_name = "nudging relaxation time scale", & units = "s", & kind = "constant scalar profile", & input_id = 1, & output_file = output_file, & grid = scalar_profile_grid, & intermediate_grid = scalar_profile_intermediate & ) output_var_table(55) % to_be_processed = ls_forcing_variables_required ! Attributes shared among all variables output_var_table(:) % source = nc_source_text END SUBROUTINE setup_variable_tables !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initializes and nc_var varible with the given parameters. The 'kind' !> parameter is used to infer the correct netCDF IDs and the level of detail, !> 'lod', as defined by the PALM-4U input data standard. !------------------------------------------------------------------------------! FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, & grid, intermediate_grid, output_file, is_profile & ) RESULT(var) CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind INTEGER, INTENT(IN) :: input_id TYPE(grid_definition), INTENT(IN), TARGET :: grid, intermediate_grid TYPE(nc_file), INTENT(IN) :: output_file LOGICAL, INTENT(IN), OPTIONAL :: is_profile CHARACTER(LEN=LNAME) :: out_var_kind TYPE(nc_var) :: var out_var_kind = TRIM(kind) IF (PRESENT(is_profile)) THEN IF (is_profile) out_var_kind = TRIM(kind) // ' profile' END IF var % name = name var % standard_name = std_name var % long_name = long_name var % units = units var % kind = TRIM(out_var_kind) var % input_id = input_id var % nt = SIZE (output_file % time) var % grid => grid var % intermediate_grid => intermediate_grid SELECT CASE( TRIM(out_var_kind) ) !TODO: Using global module variables 'init_variables_required' and !TODO: 'boundary_variables_required'. Encapsulate in settings type !TODO: and pass into init_nc_var. CASE( 'init soil' ) var % nt = 1 var % lod = 2 var % ndim = 3 var % dimids(1:3) = output_file % dimids_soil var % dimvarids(1:3) = output_file % dimvarids_soil var % to_be_processed = init_variables_required var % task = "interpolate_2d" CASE( 'init scalar' ) var % nt = 1 var % lod = 2 var % ndim = 3 var % dimids(1:3) = output_file % dimids_scl var % dimvarids(1:3) = output_file % dimvarids_scl var % to_be_processed = init_variables_required var % task = "interpolate_3d" CASE( 'init u' ) var % nt = 1 var % lod = 2 var % ndim = 3 var % dimids(1) = output_file % dimids_vel(1) var % dimids(2) = output_file % dimids_scl(2) var % dimids(3) = output_file % dimids_scl(3) var % dimvarids(1) = output_file % dimvarids_vel(1) var % dimvarids(2) = output_file % dimvarids_scl(2) var % dimvarids(3) = output_file % dimvarids_scl(3) var % to_be_processed = init_variables_required var % task = "interpolate_3d" CASE( 'init v' ) var % nt = 1 var % lod = 2 var % ndim = 3 var % dimids(1) = output_file % dimids_scl(1) var % dimids(2) = output_file % dimids_vel(2) var % dimids(3) = output_file % dimids_scl(3) var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_vel(2) var % dimvarids(3) = output_file % dimvarids_scl(3) var % to_be_processed = init_variables_required var % task = "interpolate_3d" CASE( 'init w' ) var % nt = 1 var % lod = 2 var % ndim = 3 var % dimids(1) = output_file % dimids_scl(1) var % dimids(2) = output_file % dimids_scl(2) var % dimids(3) = output_file % dimids_vel(3) var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_scl(2) var % dimvarids(3) = output_file % dimvarids_vel(3) var % to_be_processed = init_variables_required var % task = "interpolate_3d" CASE( 'init scalar profile', 'init u profile', 'init v profile') var % nt = 1 var % lod = 1 var % ndim = 1 var % dimids(1) = output_file % dimids_scl(3) !z var % dimvarids(1) = output_file % dimvarids_scl(3) !z var % to_be_processed = init_variables_required var % task = "average profile" CASE( 'init w profile') var % nt = 1 var % lod = 1 var % ndim = 1 var % dimids(1) = output_file % dimids_vel(3) !z var % dimvarids(1) = output_file % dimvarids_vel(3) !z var % to_be_processed = init_variables_required var % task = "average profile" CASE( 'surface forcing' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1:2) = output_file % dimids_soil(1:2) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1:2) = output_file % dimvarids_soil(1:2) var % to_be_processed = boundary_variables_required var % task = "interpolate_2d" CASE( 'left scalar', 'right scalar') ! same as right var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_scl(2) var % dimids(2) = output_file % dimids_scl(3) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(2) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'north scalar', 'south scalar') ! same as south var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_scl(1) var % dimids(2) = output_file % dimids_scl(3) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'top scalar', 'top w' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_scl(1) var % dimids(2) = output_file % dimids_scl(2) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_scl(2) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'left u', 'right u' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_scl(2) var % dimids(2) = output_file % dimids_scl(3) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(2) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'north u', 'south u' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_vel(1) !x var % dimids(2) = output_file % dimids_scl(3) !z var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_vel(1) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'top u' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_vel(1) !x var % dimids(2) = output_file % dimids_scl(2) !z var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_vel(1) var % dimvarids(2) = output_file % dimvarids_scl(2) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'left v', 'right v' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_vel(2) var % dimids(2) = output_file % dimids_scl(3) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_vel(2) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'north v', 'south v' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_scl(1) !x var % dimids(2) = output_file % dimids_scl(3) !z var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_scl(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'top v' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_scl(1) !x var % dimids(2) = output_file % dimids_vel(2) !z var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_vel(2) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'left w', 'right w') var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time var % dimids(1) = output_file % dimids_scl(2) var % dimids(2) = output_file % dimids_vel(3) var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(2) var % dimvarids(2) = output_file % dimvarids_vel(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'north w', 'south w' ) var % lod = -1 var % ndim = 3 var % dimids(3) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_scl(1) !x var % dimids(2) = output_file % dimids_vel(3) !z var % dimvarids(3) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(1) var % dimvarids(2) = output_file % dimvarids_vel(3) var % to_be_processed = boundary_variables_required var % task = "interpolate_3d" CASE( 'time series' ) var % lod = 0 var % ndim = 1 var % dimids(1) = output_file % dimid_time !t var % dimvarids(1) = output_file % dimvarid_time var % to_be_processed = .TRUE. var % task = "average scalar" CASE( 'constant scalar profile' ) var % lod = -1 var % ndim = 2 var % dimids(2) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_scl(3) !z var % dimvarids(2) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(3) var % to_be_processed = .TRUE. var % task = "set profile" CASE( 'large-scale scalar forcing' ) var % lod = -1 var % ndim = 2 var % dimids(2) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_scl(3) !z var % dimvarids(2) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_scl(3) var % to_be_processed = ls_forcing_variables_required var % task = "average large-scale profile" CASE( 'large-scale w forcing' ) var % lod = -1 var % ndim = 2 var % dimids(2) = output_file % dimid_time !t var % dimids(1) = output_file % dimids_vel(3) !z var % dimvarids(2) = output_file % dimvarid_time var % dimvarids(1) = output_file % dimvarids_vel(3) var % to_be_processed = ls_forcing_variables_required var % task = "average large-scale profile" CASE DEFAULT message = "Variable kind '" // TRIM(kind) // "' not recognized." CALL abort ('init_nc_var', message) END SELECT END FUNCTION init_nc_var SUBROUTINE fini_variables() CALL report('fini_variables', 'Deallocating variable table') DEALLOCATE( input_var_table ) END SUBROUTINE fini_variables SUBROUTINE fini_io_groups() CALL report('fini_io_groups', 'Deallocating IO groups') DEALLOCATE( io_group_list ) END SUBROUTINE fini_io_groups SUBROUTINE fini_file_lists() CALL report('fini_file_lists', 'Deallocating file lists') DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files ) END SUBROUTINE fini_file_lists SUBROUTINE get_input_file_list(start_date_string, start_hour, end_hour, & step_hour, path, prefix, suffix, file_list) CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, path INTEGER, INTENT(IN) :: start_hour, end_hour, step_hour CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:) INTEGER :: number_of_intervals, hour, i CHARACTER(LEN=DATE) :: date_string number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour ) ALLOCATE( file_list(number_of_intervals + 1) ) DO i = 0, number_of_intervals hour = start_hour + i * step_hour date_string = add_hours_to(start_date_string, hour) file_list(i+1) = TRIM(path) // TRIM(prefix) // TRIM(date_string) // & TRIM(suffix) // '.nc' message = "Set up input file name '" // TRIM(file_list(i+1)) //"'" CALL report('input_file_list', message) END DO END SUBROUTINE get_input_file_list !------------------------------------------------------------------------------! ! Description: ! ------------ !> Carries out any physical conversion of the quantities in the given input !> buffer needed to obtain the quantity required by PALM-4U. For instance, !> velocities are rotated to the PALM-4U coordinate system and the potential !> temperature is computed from the absolute temperature and pressure. !> !> Note, that the preprocessing does not include any grid change. The result !> array will match a COSMO-DE scalar array. !------------------------------------------------------------------------------! SUBROUTINE preprocess(group, input_buffer, cosmo_grid, iter) TYPE(io_group), INTENT(INOUT), TARGET :: group TYPE(container), INTENT(INOUT), ALLOCATABLE :: input_buffer(:) TYPE(grid_definition), INTENT(IN) :: cosmo_grid INTEGER, INTENT(IN) :: iter REAL(dp), ALLOCATABLE :: basic_state_pressure(:) TYPE(container), ALLOCATABLE :: preprocess_buffer(:) INTEGER :: hour, dt INTEGER :: i, j, k INTEGER :: nx, ny, nz input_buffer(:) % is_preprocessed = .FALSE. SELECT CASE( group % kind ) CASE( 'velocities' ) ! Allocate a compute buffer with the same number of arrays as the input ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) ) ! Allocate u and v arrays with scalar dimensions nx = SIZE(input_buffer(1) % array, 1) ny = SIZE(input_buffer(1) % array, 2) nz = SIZE(input_buffer(1) % array, 3) ALLOCATE( preprocess_buffer(1) % array(nx, ny, nz) ) ! u buffer ALLOCATE( preprocess_buffer(2) % array(nx, ny, nz) ) ! v buffer CALL run_control('time', 'alloc') ! interpolate U and V to centres CALL centre_velocities( u_face = input_buffer(1) % array, & v_face = input_buffer(2) % array, & u_centre = preprocess_buffer(1) % array, & v_centre = preprocess_buffer(2) % array ) cfg % rotation_method = 'rotated-pole' SELECT CASE(cfg % rotation_method) CASE('rotated-pole') ! rotate U and V to PALM-4U orientation and overwrite U and V with ! rotated velocities DO k = 1, nz DO j = 2, ny DO i = 2, nx CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k), & vrot = preprocess_buffer(2) % array(i,j,k), & rlat = cosmo_grid % lat(j-1), & rlon = cosmo_grid % lon(i-1), & pollat = phi_cn, & pollon = lambda_cn, & u = input_buffer(1) % array(i,j,k), & v = input_buffer(2) % array(i,j,k) ) END DO END DO END DO CASE DEFAULT message = "Rotation method '" // TRIM(cfg % rotation_method) // & "' not recognized." CALL abort('preprocess', message) END SELECT ! set values input_buffer(1) % array(1,:,:) = 0.0_dp input_buffer(2) % array(1,:,:) = 0.0_dp input_buffer(1) % array(:,1,:) = 0.0_dp input_buffer(2) % array(:,1,:) = 0.0_dp input_buffer(1:2) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') DEALLOCATE( preprocess_buffer ) CALL run_control('time', 'alloc') message = "Input buffers for group '" // TRIM(group % kind) // "'"//& " preprocessed sucessfully." CALL report('preprocess', message) CASE( 'temperature' ) ! P and T nx = SIZE(input_buffer(1) % array, 1) ny = SIZE(input_buffer(1) % array, 2) nz = SIZE(input_buffer(1) % array, 3) ! Compute absolute pressure if presure perturbation has been read in. IF ( TRIM(group % in_var_list(2) % name) == 'PP' ) THEN message = "Absolute pressure, P, not available, " // & "computing from pressure preturbation PP." CALL report('preprocess', message) ALLOCATE( basic_state_pressure(1:nz) ) CALL run_control('time', 'alloc') DO j = 1, ny DO i = 1, nx CALL get_basic_state(cosmo_grid % hfl(i,j,:), BETA, P_SL, T_SL,& RD, G, basic_state_pressure) input_buffer (2) % array(i,j,:) = & input_buffer (2) % array(i,j,:) + basic_state_pressure(:) END DO END DO CALL run_control('time', 'comp') DEALLOCATE( basic_state_pressure ) CALL run_control('time', 'alloc') group % in_var_list(2) % name = 'P' END IF ! Convert absolute to potential temperature input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) * & EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) )) input_buffer(1:2) % is_preprocessed = .TRUE. message = "Input buffers for group '" // TRIM(group % kind) // "'"//& " preprocessed sucessfully." CALL report('preprocess', message) CALL run_control('time', 'comp') CASE( 'scalar' ) ! S or W input_buffer(:) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') CASE( 'soil-temperature' ) ! CALL fill_water_cells(soiltyp, input_buffer(1) % array, & SIZE(input_buffer(1) % array, 3), & FILL_ITERATIONS) input_buffer(:) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') CASE( 'soil-water' ) ! CALL fill_water_cells(soiltyp, input_buffer(1) % array, & SIZE(input_buffer(1) % array, 3), & FILL_ITERATIONS) nx = SIZE(input_buffer(1) % array, 1) ny = SIZE(input_buffer(1) % array, 2) nz = SIZE(input_buffer(1) % array, 3) DO k = 1, nz DO j = 1, ny DO i = 1, nx input_buffer(1) % array(i,j,k) = & input_buffer(1) % array(i,j,k) * d_depth_rho_inv(k) END DO END DO END DO message = "Converted soil water from [kg/m^2] to [m^3/m^3]" CALL report('preprocess', message) input_buffer(:) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') CASE( 'surface' ) ! input_buffer(:) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') CASE( 'accumulated' ) ! message = "De-accumulating '" // TRIM(group % in_var_list(1) % name) //& "' in iteration " // TRIM(str(iter)) CALL report('preprocess', message) hour = iter - 1 dt = MODULO(hour, 3) + 1 ! averaging period SELECT CASE(dt) CASE(1) !input has been accumulated over one hour. Leave as is !input_buffer(1) % array(:,:,:) carrries one-hour integral CASE(2) !input has been accumulated over two hours. Subtract previous step !input_buffer(1) % array(:,:,:) carrries one-hour integral !input_buffer(2) % array(:,:,:) carrries two-hour integral CALL deaverage( & avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp, & avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp, & avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp ) !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour CASE(3) !input has been accumulated over three hours. Subtract previous step !input_buffer(1) % array(:,:,:) carrries three-hour integral !input_buffer(2) % array(:,:,:) still carrries two-hour integral CALL deaverage( & avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp, & avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp, & avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp ) !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA CASE DEFAULT message = "Invalid averaging period '" // TRIM(str(dt)) // " hours" CALL abort('preprocess', message) END SELECT input_buffer(:) % is_preprocessed = .TRUE. CALL run_control('time', 'comp') CASE( 'running average' ) ! message = "De-averaging '" // TRIM(group % in_var_list(1) % name) // & "' in iteration " // TRIM(str(iter)) CALL report('preprocess', message) hour = iter - 1 dt = MODULO(hour, 3) + 1 ! averaging period SELECT CASE(dt) CASE(1) !input has been accumulated over one hour. Leave as is !input_buffer(1) % array(:,:,:) carrries one-hour integral CASE(2) !input has been accumulated over two hours. Subtract previous step !input_buffer(1) % array(:,:,:) carrries one-hour integral !input_buffer(2) % array(:,:,:) carrries two-hour integral CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp, & input_buffer(2) % array(:,:,:), 2.0_dp, & input_buffer(1) % array(:,:,:), 1.0_dp) !input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour CASE(3) !input has been accumulated over three hours. Subtract previous step !input_buffer(1) % array(:,:,:) carrries three-hour integral !input_buffer(2) % array(:,:,:) still carrries two-hour integral CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp, & input_buffer(1) % array(:,:,:), 3.0_dp, & input_buffer(1) % array(:,:,:), 1.0_dp) !input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA CASE DEFAULT message = "Invalid averaging period '" // TRIM(str(dt)) // " hours" CALL abort('preprocess', message) END SELECT input_buffer(:) % is_preprocessed = .TRUE. CASE DEFAULT message = "Buffer kind '" // TRIM(group % kind) // "' is not supported." CALL abort('prerpocess', message) END SELECT CALL run_control('time', 'comp') END SUBROUTINE preprocess ! Description: ! ------------ !> Computes average soil values in COSMO-DE water cells from neighbouring !> non-water cells. This is done as a preprocessing step for the COSMO-DE !> soil input arrays, which contain unphysical values for water cells. !> !> This routine computes the average of up to all nine neighbouring cells !> or keeps the original value, if not at least one non-water neightbour !> is available. !> !> By repeatedly applying this step, soil data can be extrapolated into !> 'water' regions occupying multiple cells, one cell per iteration. !> !> Input parameters: !> ----------------- !> soiltyp : 2d map of COSMO-DE soil types !> nz : number of layers in the COSMO-DE soil !> niter : number iterations !> !> Output parameters: !> ------------------ !> array : the soil array (i.e. water content or temperature) !------------------------------------------------------------------------------! SUBROUTINE fill_water_cells(soiltyp, array, nz, niter) INTEGER(hp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: array INTEGER, INTENT(IN) :: nz, niter REAL(dp), DIMENSION(nz) :: column INTEGER(hp), DIMENSION(:,:), ALLOCATABLE :: old_soiltyp, new_soiltyp INTEGER :: l, i, j, nx, ny, n_cells, ii, jj, iter INTEGER, DIMENSION(8) :: di, dj nx = SIZE(array, 1) ny = SIZE(array, 2) di = (/ -1, -1, -1, 0, 0, 1, 1, 1 /) dj = (/ -1, 0, 1, -1, 1, -1, 0, 1 /) ALLOCATE(old_soiltyp(SIZE(soiltyp,1), & SIZE(soiltyp,2) )) ALLOCATE(new_soiltyp(SIZE(soiltyp,1), & SIZE(soiltyp,2) )) old_soiltyp(:,:) = soiltyp(:,:,1) new_soiltyp(:,:) = soiltyp(:,:,1) DO iter = 1, niter DO j = 1, ny DO i = 1, nx IF (old_soiltyp(i,j) == WATER_ID) THEN n_cells = 0 column(:) = 0.0_dp DO l = 1, SIZE(di) ii = MIN(nx, MAX(1, i + di(l))) jj = MIN(ny, MAX(1, j + dj(l))) IF (old_soiltyp(ii,jj) .NE. WATER_ID) THEN n_cells = n_cells + 1 column(:) = column(:) + array(ii,jj,:) END IF END DO ! Overwrite if at least one non-water neighbour cell is available IF (n_cells > 0) THEN array(i,j,:) = column(:) / n_cells new_soiltyp(i,j) = 0 END IF END IF END DO END DO old_soiltyp(:,:) = new_soiltyp(:,:) END DO DEALLOCATE(old_soiltyp, new_soiltyp) END SUBROUTINE fill_water_cells END MODULE grid