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

Last change on this file since 3516 was 3456, checked in by eckhard, 6 years ago

inifor: Removed surface forcing and internal arrays from netCDF output

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