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

Last change on this file since 3753 was 3680, checked in by knoop, 6 years ago

Added cpp-option netcdf to inifor

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