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

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

inifor: Re-compute geostrophic winds every time step

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