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

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

inifor: Average initial profiles over the PALM, not the geostrophic, region

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