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

Last change on this file since 3615 was 3615, checked in by raasch, 6 years ago

bugfix for last commit: abort replaced by inifor_abort

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