source: palm/trunk/UTIL/inifor/src/inifor.f90 @ 3785

Last change on this file since 3785 was 3785, checked in by eckhard, 5 years ago

inifor: Removed unused variables, improved coding style

  • Property svn:keywords set to Id
File size: 23.2 KB
RevLine 
[2696]1!> @file src/inifor.f90
2!------------------------------------------------------------------------------!
[2718]3! This file is part of the PALM model system.
[2696]4!
[2718]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
[2696]8! version.
9!
[2718]10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
[2696]13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[3779]17! Copyright 2017-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Deutscher Wetterdienst Offenbach
[2696]19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
[3183]23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor.f90 3785 2019-03-06 10:41:14Z eckhard $
[3779]28! Average geostrophic wind components on coarse COSMO levels instead of fine PALM levels
29! Remove --debug netCDF output of internal pressure profiles
30!
31! 3680 2019-01-18 14:54:12Z knoop
[3618]32! Prefixed all INIFOR modules with inifor_
33!
34!
35! 3615 2018-12-10 07:21:03Z raasch
[3615]36! bugfix: abort replaced by inifor_abort
37!
38! 3613 2018-12-07 18:20:37Z eckhard
[3613]39! Moved version output to setup_parameters()
40!
41! 3557 2018-11-22 16:01:22Z eckhard
[3557]42! Updated documentation
43!
44! 3537 2018-11-20 10:53:14Z eckhard
[3537]45! Print version number on program start
46!
47!
48! 3456 2018-10-30 14:29:54Z eckhard
[3456]49! NetCDf output of internal arrays only with --debug option
50!
51!
52! 3401 2018-10-23 09:33:31Z eckhard
[3401]53! Re-compute geostrophic winds every time step
54!
55! 3395 2018-10-22 17:32:49Z eckhard
[3395]56! Added main loop support for computation of geostrophic winds and surface
57!     pressure
58! Cleaned up terminal output, show some meVssages only with --debug option
59!
60! 3183 2018-07-27 14:25:55Z suehring
[3182]61! Introduced new PALM grid stretching
62! Renamend initial-condition mode variable 'mode' to 'ic_mode'
63! Improved log messages
[2696]64!
65!
[3183]66! 3182 2018-07-27 13:36:03Z suehring
[2696]67! Initial revision
68!
69!
70!
71! Authors:
72! --------
[3557]73!> @author Eckhard Kadasch, (Deutscher Wetterdienst, Offenbach)
[2696]74!
75! Description:
76! ------------
77!> INIFOR is an interpolation tool for generating meteorological initialization
78!> and forcing data for the urban climate model PALM-4U. The required
79!> meteorological fields are interpolated from output data of the mesoscale
[3779]80!> model COSMO. This is the main program file.
[2696]81!------------------------------------------------------------------------------!
82 PROGRAM inifor
[3680]83#if defined ( __netcdf )
[2696]84
[3618]85    USE inifor_control
86    USE inifor_defs
87    USE inifor_grid,                                                           &
[2696]88        ONLY:  setup_parameters, setup_grids, setup_variable_tables,           &
89               setup_io_groups, fini_grids, fini_variables, fini_io_groups,    &
[3182]90               fini_file_lists, preprocess, origin_lon, origin_lat,            &
[2696]91               output_file, io_group_list, output_var_table,                   &
[3395]92               cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3,                 &
93               averaging_width_ns, averaging_width_ew, phi_n, lambda_n,        &
94               lam_centre, phi_centre
[3618]95    USE inifor_io
96    USE inifor_transform,                                                      &
[3779]97        ONLY:  average_pressure_perturbation, average_profile, interpolate_1d, &
98               interpolate_1d_arr, interpolate_2d, interpolate_3d,             &
99               interp_average_profile, geostrophic_winds, extrapolate_density, &
100               extrapolate_pressure, get_surface_pressure
[3618]101    USE inifor_types
[2696]102   
103    IMPLICIT NONE
104   
[3557]105    INTEGER ::  igroup !< loop index for IO groups
106    INTEGER ::  ivar   !< loop index for output variables
107    INTEGER ::  iter   !< loop index for time steps
108
109    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
110    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
[3779]111    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
112    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
113    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
114    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
[3557]115    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
116    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
117    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
118    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
119    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
120    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
121    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
122    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
123
124    REAL(dp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
[3779]125    REAL(dp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
[3557]126
127    TYPE(nc_var), POINTER        ::  output_var      !< pointer to the currently processed output variable
128    TYPE(io_group), POINTER      ::  group           !< pointer to the currently processed IO group
129    TYPE(container), ALLOCATABLE ::  input_buffer(:) !< buffer of the current IO group's input arrays
130
131    LOGICAL, SAVE ::  ug_vg_have_been_computed = .FALSE. !< flag for managing geostrophic wind allocation and computation
[3785]132    !LOGICAL, SAVE ::  debugging_output = .TRUE.          !< flag controllging output of internal variables
[2696]133   
134!> \mainpage About INIFOR
135!>  ...
136!
137!------------------------------------------------------------------------------
138!- Section 1: Initialization
139!------------------------------------------------------------------------------
140 CALL run_control('init', 'void')
141
[3557]142!
143!-- Initialize INIFOR's parameters from command-line interface and namelists
[2696]144    CALL setup_parameters()
145
[3557]146!
147!-- Initialize all grids, including interpolation neighbours and weights
[2696]148    CALL setup_grids()
149 CALL run_control('time', 'init')
150
[3557]151!
152!-- Initialize the netCDF output file and define dimensions
[3182]153    CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
154                                 origin_lon, origin_lat)
[2696]155 CALL run_control('time', 'write')
156
[3557]157!
158!-- Set up the tables containing the input and output variables and set
159!-- the corresponding netCDF dimensions for each output variable
[3182]160    CALL setup_variable_tables(cfg % ic_mode)
[2696]161 CALL run_control('time', 'write')
162
[3557]163!
164!-- Add the output variables to the netCDF output file
[3785]165    CALL setup_netcdf_variables(output_file % name, output_var_table)
[2696]166
[3182]167    CALL setup_io_groups()
[2696]168 CALL run_control('time', 'init')
169
170!------------------------------------------------------------------------------
171!- Section 2: Main loop
172!------------------------------------------------------------------------------
173    DO igroup = 1, SIZE(io_group_list)
174
175       group => io_group_list(igroup)
176       IF ( group % to_be_processed )  THEN
177         
178          DO iter = 1, group % nt 
179
180!------------------------------------------------------------------------------
181!- Section 2.1: Read and preprocess input data
182!------------------------------------------------------------------------------
183             CALL read_input_variables(group, iter, input_buffer)
184 CALL run_control('time', 'read')
185
186             CALL preprocess(group, input_buffer, cosmo_grid, iter)
187 CALL run_control('time', 'comp')
188
[3182]189             !TODO: move this assertion into 'preprocess'.
[2696]190             IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) )  THEN
191                message = "Input buffers for group '" // TRIM(group % kind) // &
192                   "' could not be preprocessed sucessfully."
[3615]193                CALL inifor_abort('main loop', message)
[3785]194             ENDIF
[2696]195
196!------------------------------------------------------------------------------
197!- Section 2.2: Interpolate each output variable of the group
198!------------------------------------------------------------------------------
199             DO ivar = 1, group % nv
200
201                output_var => group % out_vars( ivar )
202
203                IF ( output_var % to_be_processed .AND.                        &
204                     iter .LE. output_var % nt )  THEN
205
206                   message = "Processing '" // TRIM(output_var % name) //      &
207                             "' (" // TRIM(output_var % kind) //               &
208                             "), iteration " // TRIM(str(iter)) //" of " //    &
209                             TRIM(str(output_var % nt))
210                   CALL report('main loop', message)
211
212                   SELECT CASE( TRIM(output_var % task) )
213
214                   CASE( 'interpolate_2d' ) 
215                   
216                      SELECT CASE( TRIM(output_var % kind) )
217                       
218                      CASE( 'init soil' )
219
220                         ALLOCATE( output_arr( 0:output_var % grid % nx,       &
221                                               0:output_var % grid % ny,       &
222                                               SIZE(output_var % grid % depths) ) )
223
224                      CASE ( 'surface forcing' )
225
226                         ALLOCATE( output_arr( 0:output_var % grid % nx,       &
227                                               0:output_var % grid % ny, 1 ) )
228
229                      CASE DEFAULT
230
[3182]231                          message = "'" // TRIM(output_var % kind) // "' is not a soil variable"
[3615]232                          CALL inifor_abort("main loop", message)
[2696]233
234                      END SELECT
235 CALL run_control('time', 'alloc')
236
237                      CALL interpolate_2d(input_buffer(output_var % input_id) % array(:,:,:), &
238                              output_arr(:,:,:), output_var % intermediate_grid, output_var)
239 CALL run_control('time', 'comp')
240
241
242                   CASE ( 'interpolate_3d' )
243
244                      ALLOCATE( output_arr( 0:output_var % grid % nx,          &
245                                            0:output_var % grid % ny,          &
[3182]246                                            1:output_var % grid % nz ) )
[2696]247
248 CALL run_control('time', 'alloc')
249                      CALL interpolate_3d(                                     &
250                         input_buffer(output_var % input_id) % array(:,:,:),   &
251                         output_arr(:,:,:),                                    &
252                         output_var % intermediate_grid,                       &
253                         output_var % grid)
254 CALL run_control('time', 'comp')
255
256                   CASE ( 'average profile' )
257
258                      ALLOCATE( output_arr( 0:output_var % grid % nx,          &
259                                            0:output_var % grid % ny,          &
[3182]260                                            1:output_var % grid % nz ) )
[2696]261 CALL run_control('time', 'alloc')
262                     
[3779]263                      CALL interp_average_profile(                             &
[2696]264                         input_buffer(output_var % input_id) % array(:,:,:),   &
[3395]265                         output_arr(0,0,:),                                    &
266                         output_var % averaging_grid)
267
268                      IF ( TRIM(output_var % name) ==                          &
269                           'surface_forcing_surface_pressure' )  THEN
270
271                         IF ( cfg % p0_is_set )  THEN
272                            output_arr(0,0,1) = p0
273                         ELSE
274                            CALL get_surface_pressure(                         &
275                               output_arr(0,0,:), rho_centre,                  &
276                               output_var % averaging_grid)
[3785]277                         ENDIF
[3395]278
[3785]279                      ENDIF
[2696]280 CALL run_control('time', 'comp')
281
[3395]282                   CASE ( 'internal profile' )
283
284                      message = "Averaging of internal profile for variable '" //&
285                         TRIM(output_var % name) // "' is not supported."
286
287                      SELECT CASE (TRIM(output_var % name))
288
289                      CASE('internal_density_centre')
[3779]290                         ALLOCATE( rho_centre( 1:cosmo_grid % nz) )
[3395]291                         internal_arr => rho_centre
292
293                      CASE('internal_density_north')
[3779]294                         ALLOCATE( rho_north( 1:cosmo_grid % nz) )
[3395]295                         internal_arr => rho_north
296
297                      CASE('internal_density_south')
[3779]298                         ALLOCATE( rho_south( 1:cosmo_grid % nz) )
[3395]299                         internal_arr => rho_south
300
301                      CASE('internal_density_east')
[3779]302                         ALLOCATE( rho_east( 1:cosmo_grid % nz) )
[3395]303                         internal_arr => rho_east
304
305                      CASE('internal_density_west')
[3779]306                         ALLOCATE( rho_west( 1:cosmo_grid % nz) )
[3395]307                         internal_arr => rho_west
308
309                      CASE('internal_pressure_north')
[3779]310                         ALLOCATE( p_north( 1:cosmo_grid % nz) )
[3395]311                         internal_arr => p_north
312
313                      CASE('internal_pressure_south')
[3779]314                         ALLOCATE( p_south( 1:cosmo_grid % nz) )
[3395]315                         internal_arr => p_south
316
317                      CASE('internal_pressure_east')
[3779]318                         ALLOCATE( p_east( 1:cosmo_grid % nz) )
[3395]319                         internal_arr => p_east
320
321                      CASE('internal_pressure_west')
[3779]322                         ALLOCATE( p_west( 1:cosmo_grid % nz) )
[3395]323                         internal_arr => p_west
324
325                      CASE DEFAULT
[3615]326                         CALL inifor_abort('main loop', message)
[3395]327
328                      END SELECT
329 CALL run_control('time', 'alloc')
330
331
[3779]332                      SELECT CASE( TRIM( output_var % name ) )
[3395]333
[3779]334                      CASE( 'internal_pressure_north',                         &
335                            'internal_pressure_south',                         &
336                            'internal_pressure_east',                          &
337                            'internal_pressure_west' )
[3395]338
[3779]339                         CALL average_pressure_perturbation(                   &
340                            input_buffer(output_var % input_id) % array(:,:,:),&
341                            internal_arr(:),                                   &
342                            cosmo_grid, output_var % averaging_grid            &
343                         )
[3395]344
[3779]345                      CASE DEFAULT
[3395]346
[3779]347                         CALL average_profile(                                 &
348                            input_buffer(output_var % input_id) % array(:,:,:),&
349                            internal_arr(:),                                   &
[3785]350                            output_var % averaging_grid                        &
[3779]351                         )
[3395]352
[3779]353                      END SELECT
[3395]354
355
[3779]356!
357!--                   Output of geostrophic pressure profiles (with --debug
358!--                   option) is currently deactivated, since they are now
359!--                   defined on averaged COSMO levels instead of PALM levels
360!--                   (requires definiton of COSMO levels in netCDF output.)
361                      !IF (.TRUE.)  THEN
362                      !   ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
363                      !   output_arr(1,1,:) = internal_arr(:)
[3785]364                      !ENDIF
[3395]365 CALL run_control('time', 'comp')
366
[3557]367!
368!--                This case gets called twice, the first time for ug, the
369!--                second time for vg. We compute ug and vg at the first call
[3779]370!--                and keep both of them around for the second call.
[3395]371                   CASE ( 'geostrophic winds' )
372
373                      IF (.NOT. ug_vg_have_been_computed )  THEN
[3779]374                         ALLOCATE( ug_palm(output_var % grid % nz) )
375                         ALLOCATE( vg_palm(output_var % grid % nz) )
376                         ALLOCATE( ug_cosmo(cosmo_grid % nz) )
377                         ALLOCATE( vg_cosmo(cosmo_grid % nz) )
[3395]378
[3779]379                         IF ( cfg % ug_defined_by_user )  THEN
380                            ug_palm = cfg % ug
381                            vg_palm = cfg % vg
[3395]382                         ELSE
[3557]383                            CALL geostrophic_winds( p_north, p_south, p_east,  &
384                                                    p_west, rho_centre, f3,    &
385                                                    averaging_width_ew,        &
386                                                    averaging_width_ns,        &
387                                                    phi_n, lambda_n,           &
388                                                    phi_centre, lam_centre,    &
[3779]389                                                    ug_cosmo, vg_cosmo )
390
391                            CALL interpolate_1d( ug_cosmo, ug_palm,             &
392                                                 output_var % grid )
393
394                            CALL interpolate_1d( vg_cosmo, vg_palm,             &
395                                                 output_var % grid )
[3785]396                         ENDIF
[3395]397
398                         ug_vg_have_been_computed = .TRUE.
399
[3785]400                      ENDIF
[3395]401
[3557]402!
[3779]403!--                   Select output array of current geostrophic wind component
[3395]404                      SELECT CASE(TRIM(output_var % name))
405                      CASE ('ls_forcing_ug')
[3779]406                         ug_vg_palm => ug_palm
[3395]407                      CASE ('ls_forcing_vg')
[3779]408                         ug_vg_palm => vg_palm
[3395]409                      END SELECT
410
411                      ALLOCATE( output_arr(1,1,output_var % grid % nz) )
[3779]412                      output_arr(1,1,:) = ug_vg_palm(:)
[3395]413
[2696]414                   CASE ( 'average scalar' )
415
416                      ALLOCATE( output_arr(1,1,1) )
417 CALL run_control('time', 'alloc')
418                      output_arr(1,1,1) = p0
419 CALL run_control('time', 'comp')
420
[3182]421                   CASE ( 'set profile' )
[2696]422                     
[3182]423                      ALLOCATE( output_arr( 1, 1, 1:nz ) )
[2696]424 CALL run_control('time', 'alloc')
425
426                      SELECT CASE (TRIM(output_var % name))
427
[3182]428                      CASE('nudging_tau')
429                          output_arr(1, 1, :) = NUDGING_TAU
430
[2696]431                      CASE DEFAULT
432                          message = "'" // TRIM(output_var % name) //          &
433                             "' is not a valid '" // TRIM(output_var % kind) //&
434                             "' variable kind."
[3615]435                          CALL inifor_abort('main loop', message)
[2696]436                      END SELECT
437 CALL run_control('time', 'comp')
438
[3182]439                   CASE('average large-scale profile')
440                      message = "Averaging of large-scale forcing profiles " //&
441                                "has not been implemented, yet."
[3615]442                      CALL inifor_abort('main loop', message)
[3182]443
[2696]444                   CASE DEFAULT
445                      message = "Processing task '" // TRIM(output_var % task) //&
446                               "' not recognized."
[3615]447                      CALL inifor_abort('', message)
[2696]448
449                   END SELECT
450 CALL run_control('time', 'comp')
451
452!------------------------------------------------------------------------------
453!- Section 2.3: Write current time step of current variable
454!------------------------------------------------------------------------------
[3779]455!
456!--                Output of geostrophic pressure profiles (with --debug
457!--                option) is currently deactivated, since they are now
458!--                defined on averaged COSMO levels instead of PALM levels
459!--                (requires definiton of COSMO levels in netCDF output.)
460                   !IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
461
462                   IF (.NOT. output_var % is_internal)  THEN
[3395]463                      message = "Writing variable '" // TRIM(output_var%name) // "'."
464                      CALL report('main loop', message)
[3456]465                      CALL update_output(output_var, output_arr, iter,         &
466                                         output_file, cfg)
[2696]467 CALL run_control('time', 'write')
[3785]468                   ENDIF
[2696]469
[3395]470                   IF (ALLOCATED(output_arr))  DEALLOCATE(output_arr)
[2696]471 CALL run_control('time', 'alloc')
472
[3785]473                ENDIF
[2696]474
[3557]475!
476!--          output variable loop
[3785]477             ENDDO
[2696]478
[3401]479             ug_vg_have_been_computed = .FALSE.
[3395]480             IF ( group % kind == 'thermodynamics' )  THEN
481                DEALLOCATE( rho_centre )
[3779]482                DEALLOCATE( ug_palm )
483                DEALLOCATE( vg_palm )
484                DEALLOCATE( ug_cosmo )
485                DEALLOCATE( vg_cosmo )
486                IF ( .NOT. cfg % ug_defined_by_user )  THEN
[3395]487                   DEALLOCATE( rho_north )
488                   DEALLOCATE( rho_south )
489                   DEALLOCATE( rho_east )
490                   DEALLOCATE( rho_west )
491                   DEALLOCATE( p_north )
492                   DEALLOCATE( p_south )
493                   DEALLOCATE( p_east )
494                   DEALLOCATE( p_west )
[3785]495                ENDIF
496             ENDIF
[3395]497
[3557]498!
499!--          Keep input buffer around for averaged (radiation) and
500!--          accumulated COSMO quantities (precipitation).
[2696]501             IF ( group % kind == 'running average' .OR. &
502                  group % kind == 'accumulated' )  THEN
503             ELSE
[3395]504                CALL report('main loop', 'Deallocating input buffer', cfg % debug)
[2696]505                DEALLOCATE(input_buffer)
[3785]506             ENDIF
[2696]507 CALL run_control('time', 'alloc')
508
[3557]509!
510!--       time steps / input files loop
[3785]511          ENDDO
[2696]512
513          IF (ALLOCATED(input_buffer))  THEN
[3395]514             CALL report('main loop', 'Deallocating input buffer', cfg % debug)
[2696]515             DEALLOCATE(input_buffer)
[3785]516          ENDIF
[2696]517 CALL run_control('time', 'alloc')
518
519       ELSE
520
[3182]521          message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'"
[2696]522          IF ( ALLOCATED(group % in_var_list) )  THEN
523              message = TRIM(message) // " with input variable '" //           &
524              TRIM(group % in_var_list(1) % name) // "'."
[3785]525          ENDIF
[2696]526
[3395]527          CALL report('main loop', message, cfg % debug)
[2696]528
[3557]529!
530!--    IO group % to_be_processed conditional
[3785]531       ENDIF
[2696]532
[3557]533!
534!-- IO groups loop
[3785]535    ENDDO
[2696]536
537!------------------------------------------------------------------------------
538!- Section 3: Clean up.
539!------------------------------------------------------------------------------
540    CALL fini_file_lists()
541    CALL fini_io_groups()
542    CALL fini_variables()
543    !CALL fini_grids()
544 CALL run_control('time', 'alloc')
545 CALL run_control('report', 'void')
546
[3182]547    message = "Finished writing dynamic driver '" // TRIM(output_file % name) // &
[2696]548              "' successfully."
549    CALL report('main loop', message)
550
551
[3680]552#endif
[2696]553 END PROGRAM inifor
Note: See TracBrowser for help on using the repository browser.