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

Last change on this file since 3665 was 3618, checked in by eckhard, 3 years ago

inifor: Prefixed all INIFOR modules with inifor_ and removed unused variables

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