!> @file src/inifor.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: inifor.f90 3557 2018-11-22 16:01:22Z kanani $ ! Updated documentation ! ! 3537 2018-11-20 10:53:14Z eckhard ! Print version number on program start ! ! ! 3456 2018-10-30 14:29:54Z eckhard ! NetCDf output of internal arrays only with --debug option ! ! ! 3401 2018-10-23 09:33:31Z eckhard ! Re-compute geostrophic winds every time step ! ! 3395 2018-10-22 17:32:49Z eckhard ! Added main loop support for computation of geostrophic winds and surface ! pressure ! Cleaned up terminal output, show some meVssages only with --debug option ! ! 3183 2018-07-27 14:25:55Z suehring ! Introduced new PALM grid stretching ! Renamend initial-condition mode variable 'mode' to 'ic_mode' ! Improved log messages ! ! ! 3182 2018-07-27 13:36:03Z suehring ! Initial revision ! ! ! ! Authors: ! -------- !> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach) ! ! Description: ! ------------ !> INIFOR is an interpolation tool for generating meteorological initialization !> and forcing data for the urban climate model PALM-4U. The required !> meteorological fields are interpolated from output data of the mesoscale !> model COSMO-DE. This is the main program file. !------------------------------------------------------------------------------! PROGRAM inifor USE control USE defs USE grid, & ONLY: setup_parameters, setup_grids, setup_variable_tables, & setup_io_groups, fini_grids, fini_variables, fini_io_groups, & fini_file_lists, preprocess, origin_lon, origin_lat, & output_file, io_group_list, output_var_table, & cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3, & averaging_width_ns, averaging_width_ew, phi_n, lambda_n, & lam_centre, phi_centre USE io USE transform, & ONLY: average_profile, interpolate_2d, interpolate_3d, & geostrophic_winds, extrapolate_density, extrapolate_pressure, & get_surface_pressure USE types IMPLICIT NONE INTEGER :: igroup !< loop index for IO groups INTEGER :: ivar !< loop index for output variables INTEGER :: iter !< loop index for time steps REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr !< array buffer for interpolated quantities REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre !< density profile of the centre averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_arr !< geostrophic wind in x direction REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_arr !< geostrophic wind in y direction REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north !< density profile of the northern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south !< density profile of the southern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east !< density profile of the eastern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west !< density profile of the western averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north !< pressure profile of the northern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south !< pressure profile of the southern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east !< pressure profile of the eastern averaging domain REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west !< pressure profile of the western averaging domain REAL(dp), POINTER, DIMENSION(:) :: internal_arr !< pointer to the currently processed internal array (density, pressure) REAL(dp), POINTER, DIMENSION(:) :: ug_vg_arr !< pointer to the currently processed geostrophic wind component TYPE(nc_var), POINTER :: output_var !< pointer to the currently processed output variable TYPE(io_group), POINTER :: group !< pointer to the currently processed IO group TYPE(container), ALLOCATABLE :: input_buffer(:) !< buffer of the current IO group's input arrays LOGICAL, SAVE :: ug_vg_have_been_computed = .FALSE. !< flag for managing geostrophic wind allocation and computation LOGICAL, SAVE :: debugging_output = .TRUE. !< flag controllging output of internal variables !> \mainpage About INIFOR !> ... ! !------------------------------------------------------------------------------ !- Section 1: Initialization !------------------------------------------------------------------------------ CALL run_control('init', 'void') CALL report('main_loop', 'Running INIFOR version ' // VERSION) ! !-- Initialize INIFOR's parameters from command-line interface and namelists CALL setup_parameters() ! !-- Initialize all grids, including interpolation neighbours and weights CALL setup_grids() CALL run_control('time', 'init') ! !-- Initialize the netCDF output file and define dimensions CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date, & origin_lon, origin_lat) CALL run_control('time', 'write') ! !-- Set up the tables containing the input and output variables and set !-- the corresponding netCDF dimensions for each output variable CALL setup_variable_tables(cfg % ic_mode) CALL run_control('time', 'write') ! !-- Add the output variables to the netCDF output file CALL setup_netcdf_variables(output_file % name, output_var_table, & cfg % debug) CALL setup_io_groups() CALL run_control('time', 'init') !------------------------------------------------------------------------------ !- Section 2: Main loop !------------------------------------------------------------------------------ DO igroup = 1, SIZE(io_group_list) group => io_group_list(igroup) IF ( group % to_be_processed ) THEN DO iter = 1, group % nt !------------------------------------------------------------------------------ !- Section 2.1: Read and preprocess input data !------------------------------------------------------------------------------ CALL read_input_variables(group, iter, input_buffer) CALL run_control('time', 'read') CALL preprocess(group, input_buffer, cosmo_grid, iter) CALL run_control('time', 'comp') !TODO: move this assertion into 'preprocess'. IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) ) THEN message = "Input buffers for group '" // TRIM(group % kind) // & "' could not be preprocessed sucessfully." CALL abort('main loop', message) END IF !------------------------------------------------------------------------------ !- Section 2.2: Interpolate each output variable of the group !------------------------------------------------------------------------------ DO ivar = 1, group % nv output_var => group % out_vars( ivar ) IF ( output_var % to_be_processed .AND. & iter .LE. output_var % nt ) THEN message = "Processing '" // TRIM(output_var % name) // & "' (" // TRIM(output_var % kind) // & "), iteration " // TRIM(str(iter)) //" of " // & TRIM(str(output_var % nt)) CALL report('main loop', message) SELECT CASE( TRIM(output_var % task) ) CASE( 'interpolate_2d' ) SELECT CASE( TRIM(output_var % kind) ) CASE( 'init soil' ) ALLOCATE( output_arr( 0:output_var % grid % nx, & 0:output_var % grid % ny, & SIZE(output_var % grid % depths) ) ) CASE ( 'surface forcing' ) ALLOCATE( output_arr( 0:output_var % grid % nx, & 0:output_var % grid % ny, 1 ) ) CASE DEFAULT message = "'" // TRIM(output_var % kind) // "' is not a soil variable" CALL abort("main loop", message) END SELECT CALL run_control('time', 'alloc') CALL interpolate_2d(input_buffer(output_var % input_id) % array(:,:,:), & output_arr(:,:,:), output_var % intermediate_grid, output_var) CALL run_control('time', 'comp') CASE ( 'interpolate_3d' ) ALLOCATE( output_arr( 0:output_var % grid % nx, & 0:output_var % grid % ny, & 1:output_var % grid % nz ) ) CALL run_control('time', 'alloc') CALL interpolate_3d( & input_buffer(output_var % input_id) % array(:,:,:), & output_arr(:,:,:), & output_var % intermediate_grid, & output_var % grid) CALL run_control('time', 'comp') CASE ( 'average profile' ) ALLOCATE( output_arr( 0:output_var % grid % nx, & 0:output_var % grid % ny, & 1:output_var % grid % nz ) ) CALL run_control('time', 'alloc') CALL average_profile( & input_buffer(output_var % input_id) % array(:,:,:), & output_arr(0,0,:), & output_var % averaging_grid) IF ( TRIM(output_var % name) == & 'surface_forcing_surface_pressure' ) THEN IF ( cfg % p0_is_set ) THEN output_arr(0,0,1) = p0 ELSE CALL get_surface_pressure( & output_arr(0,0,:), rho_centre, & output_var % averaging_grid) END IF END IF CALL run_control('time', 'comp') CASE ( 'internal profile' ) message = "Averaging of internal profile for variable '" //& TRIM(output_var % name) // "' is not supported." SELECT CASE (TRIM(output_var % name)) CASE('internal_density_centre') ALLOCATE( rho_centre( 1:output_var % grid % nz ) ) internal_arr => rho_centre CASE('internal_density_north') ALLOCATE( rho_north( 1:output_var % grid % nz ) ) internal_arr => rho_north CASE('internal_density_south') ALLOCATE( rho_south( 1:output_var % grid % nz ) ) internal_arr => rho_south CASE('internal_density_east') ALLOCATE( rho_east( 1:output_var % grid % nz) ) internal_arr => rho_east CASE('internal_density_west') ALLOCATE( rho_west( 1:output_var % grid % nz ) ) internal_arr => rho_west CASE('internal_pressure_north') ALLOCATE( p_north( 1:output_var % grid % nz ) ) internal_arr => p_north CASE('internal_pressure_south') ALLOCATE( p_south( 1:output_var % grid % nz ) ) internal_arr => p_south CASE('internal_pressure_east') ALLOCATE( p_east( 1:output_var % grid % nz) ) internal_arr => p_east CASE('internal_pressure_west') ALLOCATE( p_west( 1:output_var % grid % nz ) ) internal_arr => p_west CASE DEFAULT CALL abort('main loop', message) END SELECT CALL run_control('time', 'alloc') CALL average_profile( & input_buffer(output_var % input_id) % array(:,:,:),& internal_arr(:), & output_var % averaging_grid) SELECT CASE (TRIM(output_var % name)) CASE('internal_density_centre', & 'internal_density_north', & 'internal_density_south', & 'internal_density_east', & 'internal_density_west') CALL extrapolate_density(internal_arr, & output_var % averaging_grid) CASE('internal_pressure_north') CALL extrapolate_pressure(internal_arr, rho_north, & output_var % averaging_grid) CASE('internal_pressure_south') CALL extrapolate_pressure(internal_arr, rho_south, & output_var % averaging_grid) CASE('internal_pressure_east') CALL extrapolate_pressure(internal_arr, rho_east, & output_var % averaging_grid) CASE('internal_pressure_west') CALL extrapolate_pressure(internal_arr, rho_west, & output_var % averaging_grid) CASE DEFAULT CALL abort('main loop', message) END SELECT IF (.TRUE.) THEN ALLOCATE( output_arr(1,1,1:output_var % grid % nz) ) output_arr(1,1,:) = internal_arr(:) END IF CALL run_control('time', 'comp') ! !-- This case gets called twice, the first time for ug, the !-- second time for vg. We compute ug and vg at the first call !-- and keep vg (and ug for that matter) around for the second !-- call. CASE ( 'geostrophic winds' ) IF (.NOT. ug_vg_have_been_computed ) THEN ALLOCATE( ug_arr(output_var % grid % nz) ) ALLOCATE( vg_arr(output_var % grid % nz) ) IF ( cfg % ug_is_set ) THEN ug_arr = cfg % ug vg_arr = cfg % vg ELSE CALL geostrophic_winds( p_north, p_south, p_east, & p_west, rho_centre, f3, & averaging_width_ew, & averaging_width_ns, & phi_n, lambda_n, & phi_centre, lam_centre, & ug_arr, vg_arr ) END IF ug_vg_have_been_computed = .TRUE. END IF ! !-- Prepare output of geostrophic winds SELECT CASE(TRIM(output_var % name)) CASE ('ls_forcing_ug') ug_vg_arr => ug_arr CASE ('ls_forcing_vg') ug_vg_arr => vg_arr END SELECT ALLOCATE( output_arr(1,1,output_var % grid % nz) ) output_arr(1,1,:) = ug_vg_arr(:) CASE ( 'average scalar' ) ALLOCATE( output_arr(1,1,1) ) CALL run_control('time', 'alloc') output_arr(1,1,1) = p0 CALL run_control('time', 'comp') CASE ( 'set profile' ) ALLOCATE( output_arr( 1, 1, 1:nz ) ) CALL run_control('time', 'alloc') SELECT CASE (TRIM(output_var % name)) CASE('nudging_tau') output_arr(1, 1, :) = NUDGING_TAU CASE DEFAULT message = "'" // TRIM(output_var % name) // & "' is not a valid '" // TRIM(output_var % kind) //& "' variable kind." CALL abort('main loop', message) END SELECT CALL run_control('time', 'comp') CASE('average large-scale profile') message = "Averaging of large-scale forcing profiles " //& "has not been implemented, yet." CALL abort('main loop', message) CASE DEFAULT message = "Processing task '" // TRIM(output_var % task) //& "' not recognized." CALL abort('', message) END SELECT CALL run_control('time', 'comp') !------------------------------------------------------------------------------ !- Section 2.3: Write current time step of current variable !------------------------------------------------------------------------------ IF (.NOT. output_var % is_internal .OR. debugging_output) THEN message = "Writing variable '" // TRIM(output_var%name) // "'." CALL report('main loop', message) CALL update_output(output_var, output_arr, iter, & output_file, cfg) CALL run_control('time', 'write') END IF IF (ALLOCATED(output_arr)) DEALLOCATE(output_arr) CALL run_control('time', 'alloc') END IF ! !-- output variable loop END DO ug_vg_have_been_computed = .FALSE. IF ( group % kind == 'thermodynamics' ) THEN DEALLOCATE( rho_centre ) DEALLOCATE( ug_arr ) DEALLOCATE( vg_arr ) IF ( .NOT. cfg % ug_is_set ) THEN DEALLOCATE( rho_north ) DEALLOCATE( rho_south ) DEALLOCATE( rho_east ) DEALLOCATE( rho_west ) DEALLOCATE( p_north ) DEALLOCATE( p_south ) DEALLOCATE( p_east ) DEALLOCATE( p_west ) END IF END IF ! !-- Keep input buffer around for averaged (radiation) and !-- accumulated COSMO quantities (precipitation). IF ( group % kind == 'running average' .OR. & group % kind == 'accumulated' ) THEN ELSE CALL report('main loop', 'Deallocating input buffer', cfg % debug) DEALLOCATE(input_buffer) END IF CALL run_control('time', 'alloc') ! !-- time steps / input files loop END DO IF (ALLOCATED(input_buffer)) THEN CALL report('main loop', 'Deallocating input buffer', cfg % debug) DEALLOCATE(input_buffer) END IF CALL run_control('time', 'alloc') ELSE message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'" IF ( ALLOCATED(group % in_var_list) ) THEN message = TRIM(message) // " with input variable '" // & TRIM(group % in_var_list(1) % name) // "'." END IF CALL report('main loop', message, cfg % debug) ! !-- IO group % to_be_processed conditional END IF ! !-- IO groups loop END DO !------------------------------------------------------------------------------ !- Section 3: Clean up. !------------------------------------------------------------------------------ CALL fini_file_lists() CALL fini_io_groups() CALL fini_variables() !CALL fini_grids() CALL run_control('time', 'alloc') CALL run_control('report', 'void') message = "Finished writing dynamic driver '" // TRIM(output_file % name) // & "' successfully." CALL report('main loop', message) END PROGRAM inifor