!> @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-2020 Leibniz Universitaet Hannover ! Copyright 2017-2020 Deutscher Wetterdienst Offenbach !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: inifor.f90 4523 2020-05-07 15:58:16Z eckhard $ ! respect integer working precision (iwp) specified in inifor_defs.f90 ! ! ! 4481 2020-03-31 18:55:54Z maronga ! Alert user to warnings if any ! ! ! 3866 2019-04-05 14:25:01Z eckhard ! Use PALM's working precision ! Show error message if compiled without netCDF support ! Renamed run_control -> log_runtime ! Improved coding style added comments ! ! ! 3785 2019-03-06 10:41:14Z eckhard ! Average geostrophic wind components on coarse COSMO levels instead of fine PALM levels ! Remove --debug netCDF output of internal pressure profiles ! ! 3680 2019-01-18 14:54:12Z knoop ! Prefixed all INIFOR modules with inifor_ ! ! ! 3615 2018-12-10 07:21:03Z raasch ! bugfix: abort replaced by inifor_abort ! ! 3613 2018-12-07 18:20:37Z eckhard ! Moved version output to setup_parameters() ! ! 3557 2018-11-22 16:01:22Z eckhard ! 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. This is the main program file. !------------------------------------------------------------------------------! PROGRAM inifor #if defined ( __netcdf ) USE inifor_control USE inifor_defs USE inifor_grid, & ONLY: averaging_width_ns, & averaging_width_ew, & cfg, & cosmo_grid, & f3, & fini_grids, & fini_io_groups, & fini_variables, & fini_file_lists, & io_group_list, & lam_centre, & lambda_n, & nx, ny, nz, & origin_lat, & origin_lon, & output_file, & output_var_table, & p0, & phi_centre, & phi_n, & preprocess, & palm_grid, & setup_grids, & setup_parameters, & setup_variable_tables, & setup_io_groups USE inifor_io USE inifor_transform, & ONLY: average_pressure_perturbation, & average_profile, & geostrophic_winds, & get_surface_pressure, & interp_average_profile, & interpolate_1d, & interpolate_1d_arr, & interpolate_2d, & interpolate_3d USE inifor_types IMPLICIT NONE INTEGER(iwp) :: igroup !< loop index for IO groups INTEGER(iwp) :: ivar !< loop index for output variables INTEGER(iwp) :: iter !< loop index for time steps REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: output_arr !< array buffer for interpolated quantities REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_centre !< density profile of the centre averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_cosmo !< profile of the geostrophic wind in x direction on COSMO levels REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_cosmo !< profile of the geostrophic wind in y direction on COSMO levels REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: ug_palm !< profile of the geostrophic wind in x direction interpolated onto PALM levels REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: vg_palm !< profile of the geostrophic wind in y direction interpolated onto PALM levels REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_north !< density profile of the northern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_south !< density profile of the southern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_east !< density profile of the eastern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: rho_west !< density profile of the western averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: p_north !< pressure profile of the northern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: p_south !< pressure profile of the southern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: p_east !< pressure profile of the eastern averaging domain REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET :: p_west !< pressure profile of the western averaging domain REAL(wp), POINTER, DIMENSION(:) :: internal_arr !< pointer to the currently processed internal array (density, pressure) REAL(wp), POINTER, DIMENSION(:) :: ug_vg_palm !< 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 log_runtime( 'init', 'void' ) ! !-- 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 log_runtime( '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 log_runtime( '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 log_runtime( 'time', 'write' ) ! !-- Add the output variables to the netCDF output file CALL setup_netcdf_variables( output_file%name, output_var_table ) CALL setup_io_groups CALL log_runtime( 'time', 'init' ) !------------------------------------------------------------------------------ !-- Section 2: Main loop !------------------------------------------------------------------------------ ! !-- Input and output variables are organized into IO groups. For instance, the !-- 'thermodynamics' IO group bundles the input variaebls T, P, QV and the !-- output variables p, theta, rho, and qv. !-- An IO group bunldes variables that are physically dependent on each other. !-- In case of the 'thermodynamics' group, theta = f(P,T), rho = f(P,T,QV). DO igroup = 1, SIZE( io_group_list ) group => io_group_list(igroup) IF ( group%to_be_processed ) THEN !-- Loop over all output time steps for the current group. DO iter = 1, group%nt !------------------------------------------------------------------------------ !-- Section 2.1: Read and preprocess input data !------------------------------------------------------------------------------ CALL read_input_variables( group, iter, input_buffer ) CALL log_runtime( 'time', 'read' ) !-- Carry out all required physical conversion of the input variables !-- of the current IO group on the input (COSMO) grid. For instance, !-- horizontal velocities are rotated to the PALM coordinate system and !-- potential temperature is computed from the absolute temperature and !-- pressure. CALL preprocess( group, input_buffer, cosmo_grid, iter ) CALL log_runtime( '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 inifor_abort( 'main loop', message ) ENDIF !------------------------------------------------------------------------------ !-- 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 ) ) !-- 2D horizontal interpolation 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 inifor_abort( "main loop", message ) END SELECT CALL log_runtime( 'time', 'alloc' ) CALL interpolate_2d( input_buffer(output_var%input_id)%array(:,:,:), & output_arr(:,:,:), output_var%intermediate_grid, output_var ) CALL log_runtime( 'time', 'comp' ) !-- Interpolation in 3D, used for atmospheric initial and !-- boundary conditions. CASE ( 'interpolate_3d' ) ALLOCATE( output_arr(0:output_var%grid % nx, & 0:output_var%grid % ny, & 1:output_var%grid % nz) ) CALL log_runtime( 'time', 'alloc' ) CALL interpolate_3d( & input_buffer(output_var%input_id)%array(:,:,:), & output_arr(:,:,:), & output_var%intermediate_grid, & output_var%grid) CALL log_runtime( 'time', 'comp' ) !-- Compute initial avaerage profiles (if --init-mode profile !-- is used) CASE ( 'average profile' ) ALLOCATE( output_arr(0:output_var%grid%nx, & 0:output_var%grid%ny, & 1:output_var%grid%nz) ) CALL log_runtime( 'time', 'alloc' ) CALL interp_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 ) ENDIF ENDIF CALL log_runtime( 'time', 'comp' ) !-- Compute internal profiles, required for differentiation of !-- geostrophic wind 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:cosmo_grid%nz) ) internal_arr => rho_centre CASE( 'internal_density_north' ) ALLOCATE( rho_north(1:cosmo_grid%nz) ) internal_arr => rho_north CASE( 'internal_density_south' ) ALLOCATE( rho_south(1:cosmo_grid%nz) ) internal_arr => rho_south CASE( 'internal_density_east' ) ALLOCATE( rho_east(1:cosmo_grid%nz) ) internal_arr => rho_east CASE( 'internal_density_west' ) ALLOCATE( rho_west(1:cosmo_grid%nz) ) internal_arr => rho_west CASE( 'internal_pressure_north' ) ALLOCATE( p_north(1:cosmo_grid%nz) ) internal_arr => p_north CASE( 'internal_pressure_south' ) ALLOCATE( p_south(1:cosmo_grid%nz) ) internal_arr => p_south CASE( 'internal_pressure_east' ) ALLOCATE( p_east(1:cosmo_grid%nz) ) internal_arr => p_east CASE( 'internal_pressure_west' ) ALLOCATE( p_west(1:cosmo_grid%nz) ) internal_arr => p_west CASE DEFAULT CALL inifor_abort( 'main loop', message ) END SELECT CALL log_runtime( 'time', 'alloc' ) SELECT CASE( TRIM( output_var%name ) ) CASE( 'internal_pressure_north', & 'internal_pressure_south', & 'internal_pressure_east', & 'internal_pressure_west' ) CALL average_pressure_perturbation( & input_buffer(output_var%input_id) % array(:,:,:),& internal_arr(:), & cosmo_grid, output_var%averaging_grid & ) CASE DEFAULT CALL average_profile( & input_buffer(output_var%input_id) % array(:,:,:),& internal_arr(:), & output_var%averaging_grid & ) END SELECT ! !-- Output of geostrophic pressure profiles (with --debug !-- option) is currently deactivated, since they are now !-- defined on averaged COSMO levels instead of PALM levels !-- (requires definiton of COSMO levels in netCDF output.) !IF (.TRUE.) THEN ! ALLOCATE( output_arr(1,1,1:output_var%grid % nz) ) ! output_arr(1,1,:) = internal_arr(:) !ENDIF CALL log_runtime( '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 both of them around for the second call. CASE ( 'geostrophic winds' ) IF (.NOT. ug_vg_have_been_computed ) THEN ALLOCATE( ug_palm(output_var%grid%nz) ) ALLOCATE( vg_palm(output_var%grid%nz) ) ALLOCATE( ug_cosmo(cosmo_grid%nz) ) ALLOCATE( vg_cosmo(cosmo_grid%nz) ) IF ( cfg%ug_defined_by_user ) THEN ug_palm = cfg%ug vg_palm = 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_cosmo, vg_cosmo ) CALL interpolate_1d( ug_cosmo, ug_palm, & output_var%grid ) CALL interpolate_1d( vg_cosmo, vg_palm, & output_var%grid ) ENDIF ug_vg_have_been_computed = .TRUE. ENDIF ! !-- Select output array of current geostrophic wind component SELECT CASE( TRIM( output_var%name ) ) CASE ( 'ls_forcing_ug' ) ug_vg_palm => ug_palm CASE ( 'ls_forcing_vg' ) ug_vg_palm => vg_palm END SELECT ALLOCATE( output_arr(1,1,output_var%grid%nz) ) output_arr(1,1,:) = ug_vg_palm(:) !-- User defined constant profiles CASE ( 'set profile' ) ALLOCATE( output_arr(1,1,1:nz) ) CALL log_runtime( '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 inifor_abort( 'main loop', message ) END SELECT CALL log_runtime( 'time', 'comp' ) CASE DEFAULT message = "Processing task '" // TRIM( output_var%task ) //& "' not recognized." CALL inifor_abort( '', message ) END SELECT CALL log_runtime( 'time', 'comp' ) !------------------------------------------------------------------------------ !- Section 2.3: Write current time step of current variable !------------------------------------------------------------------------------ ! !-- Output of geostrophic pressure profiles (with --debug !-- option) is currently deactivated, since they are now !-- defined on averaged COSMO levels instead of PALM levels !-- (requires definiton of COSMO levels in netCDF output.) !IF (.NOT. output_var%is_internal .OR. debugging_output) THEN IF ( .NOT. output_var%is_internal ) 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 log_runtime( 'time', 'write' ) ENDIF IF ( ALLOCATED( output_arr ) ) DEALLOCATE( output_arr ) CALL log_runtime( 'time', 'alloc' ) ENDIF ! !-- output variable loop ENDDO ug_vg_have_been_computed = .FALSE. IF ( group%kind == 'thermodynamics' ) THEN DEALLOCATE( rho_centre ) DEALLOCATE( ug_palm ) DEALLOCATE( vg_palm ) DEALLOCATE( ug_cosmo ) DEALLOCATE( vg_cosmo ) IF ( .NOT. cfg%ug_defined_by_user ) 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 ) ENDIF ENDIF ! !-- 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 ) ENDIF CALL log_runtime( 'time', 'alloc' ) ! !-- time steps / input files loop ENDDO IF ( ALLOCATED( input_buffer ) ) THEN CALL report( 'main loop', 'Deallocating input buffer', cfg%debug ) DEALLOCATE( input_buffer ) ENDIF CALL log_runtime( '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 ) // "'." ENDIF CALL report( 'main loop', message, cfg%debug ) ! !-- IO group%to_be_processed conditional ENDIF ! !-- IO groups loop ENDDO !------------------------------------------------------------------------------ !- Section 3: Clean up. !------------------------------------------------------------------------------ CALL fini_file_lists CALL fini_io_groups CALL fini_variables !CALL fini_grids CALL log_runtime( 'time', 'alloc' ) CALL log_runtime( 'report', 'void' ) message = "Finished writing dynamic driver '" // TRIM( output_file%name ) IF ( n_wrngs > 0 ) THEN message = TRIM( message ) // "' with " // TRIM( str( n_wrngs ) ) // & " warning(s). Please see logfile '" // LOG_FILE_NAME // "'." ELSE message = TRIM( message ) // "' successfully." END IF CALL report( 'main loop', message ) CALL close_log #else USE inifor_control IMPLICIT NONE message = "INIFOR was compiled without netCDF support, which is required for it to run. " // & "To use INIFOR, recompile PALM with netCDF support by adding the -D__netcdf " // & "precompiler flag to your .palm.config file." CALL inifor_abort( 'main loop', message ) #endif END PROGRAM inifor