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

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

inifor: Added computation of geostrophic winds from COSMO input

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