source: palm/trunk/SOURCE/virtual_flight_mod.f90

Last change on this file was 4843, checked in by raasch, 5 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

  • Property svn:keywords set to Id
File size: 46.1 KB
Line 
1!> @file virtual_flights_mod.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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: virtual_flight_mod.f90 4843 2021-01-15 15:22:11Z banzhafs $
27! local namelist parameter added to switch off the module although the respective module namelist
28! appears in the namelist file
29!
30! 4842 2021-01-14 10:42:28Z raasch
31! reading of namelist file and actions in case of namelist errors revised so that statement labels
32! and goto statements are not required any more
33!
34! 4828 2021-01-05 11:21:41Z Giersch
35! Bugfix, use time_since_reference_point instead of simulated_time
36!
37! 4535 2020-05-15 12:07:23Z raasch
38! bugfix for restart data format query
39!
40! 4522 2020-05-06 14:17:05Z suehring
41! Modularize user_init_flight in order to provide an explicit interface.
42!
43! 4497 2020-04-15 10:20:51Z raasch
44! file re-formatted to follow the PALM coding standard
45!
46!
47! 4495 2020-04-13 20:11:20Z raasch
48! restart data handling with MPI-IO added
49!
50! 4360 2020-01-07 11:25:50Z suehring
51! Corrected "Former revisions" section
52!
53! 4004 2019-05-24 11:32:38Z suehring
54! Allow variable start- and end locations also in return mode
55!
56! 3885 2019-04-11 11:29:34Z kanani
57! Changes related to global restructuring of location messages and introduction of additional
58! debug messages
59!
60! 3655 2019-01-07 16:51:22Z knoop
61! variables documented
62!
63! 1957 2016-07-07 10:43:48Z suehring
64! Initial revision
65!
66! Description:
67! ------------
68!> Module for virtual flight measurements.
69!> @todo Err msg PA0438: flight can be inside topography -> extra check?
70!--------------------------------------------------------------------------------------------------!
71 MODULE flight_mod
72
73    USE control_parameters,                                                                        &
74        ONLY:  debug_output,                                                                       &
75               fl_max, num_leg,                                                                    &
76               num_var_fl,                                                                         &
77               num_var_fl_user,                                                                    &
78               restart_data_format_output,                                                         &
79               virtual_flight
80
81    USE kinds
82
83    USE restart_data_mpi_io_mod,                                                                   &
84        ONLY:  rd_mpi_io_check_array, rrd_mpi_io_global_array, wrd_mpi_io_global_array
85
86    USE  user_init_flight_mod,                                                                     &
87        ONLY:  user_init_flight
88
89
90    CHARACTER(LEN=6), DIMENSION(fl_max) ::  leg_mode = 'cyclic'  !< flight mode through the model domain, either 'cyclic' or
91                                                                 !<'return'
92
93    INTEGER(iwp) ::  l           !< index for flight leg
94    INTEGER(iwp) ::  var_index   !< index for measured variable
95
96    LOGICAL, DIMENSION(:), ALLOCATABLE  ::  cyclic_leg  !< flag to identify fly mode
97
98    REAL(wp) ::  flight_end = 9999999.9_wp  !< end time of virtual flight
99    REAL(wp) ::  flight_begin = 0.0_wp      !< end time of virtual flight
100
101    REAL(wp), DIMENSION(fl_max) ::  flight_angle = 45.0_wp    !< angle determining the horizontal flight direction
102    REAL(wp), DIMENSION(fl_max) ::  flight_level = 100.0_wp   !< flight level
103    REAL(wp), DIMENSION(fl_max) ::  max_elev_change = 0.0_wp  !< maximum elevation change for the respective flight leg
104    REAL(wp), DIMENSION(fl_max) ::  rate_of_climb = 0.0_wp    !< rate of climb or descent
105    REAL(wp), DIMENSION(fl_max) ::  speed_agl = 25.0_wp       !< absolute horizontal flight speed above ground level (agl)
106    REAL(wp), DIMENSION(fl_max) ::  x_start = 999999999.0_wp  !< start x position
107    REAL(wp), DIMENSION(fl_max) ::  x_end   = 999999999.0_wp  !< end x position
108    REAL(wp), DIMENSION(fl_max) ::  y_start = 999999999.0_wp  !< start y position
109    REAL(wp), DIMENSION(fl_max) ::  y_end   = 999999999.0_wp  !< end y position
110
111    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_agl  !< u-component of flight speed
112    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_agl  !< v-component of flight speed
113    REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_agl  !< w-component of flight speed
114    REAL(wp), DIMENSION(:), ALLOCATABLE ::  x_pos  !< aircraft x-position
115    REAL(wp), DIMENSION(:), ALLOCATABLE ::  y_pos  !< aircraft y-position
116    REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_pos  !< aircraft z-position
117
118    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor_l  !< measured data on local PE
119    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor    !< measured data
120
121    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var_u  !< dummy array for possibly user-defined quantities
122
123    SAVE
124
125    PRIVATE
126
127    INTERFACE flight_header
128       MODULE PROCEDURE flight_header
129    END INTERFACE flight_header
130
131    INTERFACE flight_init
132       MODULE PROCEDURE flight_init
133    END INTERFACE flight_init
134
135    INTERFACE flight_init_output
136       MODULE PROCEDURE flight_init_output
137    END INTERFACE flight_init_output
138
139    INTERFACE flight_check_parameters
140       MODULE PROCEDURE flight_check_parameters
141    END INTERFACE flight_check_parameters
142
143    INTERFACE flight_parin
144       MODULE PROCEDURE flight_parin
145    END INTERFACE flight_parin
146
147    INTERFACE interpolate_xyz
148       MODULE PROCEDURE interpolate_xyz
149    END INTERFACE interpolate_xyz
150
151    INTERFACE flight_measurement
152       MODULE PROCEDURE flight_measurement
153    END INTERFACE flight_measurement
154
155    INTERFACE flight_rrd_global
156       MODULE PROCEDURE flight_rrd_global_ftn
157       MODULE PROCEDURE flight_rrd_global_mpi
158    END INTERFACE flight_rrd_global
159
160    INTERFACE flight_wrd_global
161       MODULE PROCEDURE flight_wrd_global
162    END INTERFACE flight_wrd_global
163
164!
165!-- Private interfaces
166    PRIVATE flight_check_parameters,                                                               &
167            flight_init_output,                                                                    &
168            interpolate_xyz
169!
170!-- Public interfaces
171    PUBLIC flight_init,                                                                            &
172           flight_header,                                                                          &
173           flight_parin,                                                                           &
174           flight_measurement,                                                                     &
175           flight_wrd_global,                                                                      &
176           flight_rrd_global
177!
178!-- Public variables
179    PUBLIC fl_max,                                                                                 &
180           sensor,                                                                                 &
181           x_pos,                                                                                  &
182           y_pos,                                                                                  &
183           z_pos
184
185 CONTAINS
186
187!--------------------------------------------------------------------------------------------------!
188! Description:
189! ------------
190!> Header output for flight module.
191!--------------------------------------------------------------------------------------------------!
192 SUBROUTINE flight_header ( io )
193
194
195    IMPLICIT NONE
196
197    INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
198
199    WRITE ( io, 1  )
200    WRITE ( io, 2  )
201    WRITE ( io, 3  ) num_leg
202    WRITE ( io, 4  ) flight_begin
203    WRITE ( io, 5  ) flight_end
204
205    DO  l=1, num_leg
206       WRITE ( io, 6   ) l
207       WRITE ( io, 7   ) speed_agl(l)
208       WRITE ( io, 8   ) flight_level(l)
209       WRITE ( io, 9   ) max_elev_change(l)
210       WRITE ( io, 10  ) rate_of_climb(l)
211       WRITE ( io, 11  ) leg_mode(l)
212    ENDDO
213
214
215 1   FORMAT (' Virtual flights: ----------------' )
216 2   FORMAT ('       Output every timestep' )
217 3   FORMAT ('       Number of flight legs:',    I3 )
218 4   FORMAT ('       Begin of measurements:',    F10.1    , ' s' )
219 5   FORMAT ('       End of measurements:',      F10.1    , ' s' )
220 6   FORMAT ('       Leg', I3/, '       ------' )
221 7   FORMAT ('          Flight speed            : ', F5.1, ' m/s' )
222 8   FORMAT ('          Flight level            : ', F5.1, ' m' )
223 9   FORMAT ('          Maximum elevation change: ', F5.1, ' m/s' )
224 10  FORMAT ('          Rate of climb / descent : ', F5.1, ' m/s' )
225 11  FORMAT ('          Leg mode                : ', A/ )
226
227 END SUBROUTINE flight_header
228
229!--------------------------------------------------------------------------------------------------!
230! Description:
231! ------------
232!> Reads the namelist flight_par.
233!--------------------------------------------------------------------------------------------------!
234 SUBROUTINE flight_parin
235
236    IMPLICIT NONE
237
238    CHARACTER(LEN=100) ::  line  !< dummy string that contains the current line of the parameter file
239
240    INTEGER(iwp) ::  io_status   !< status after reading the namelist file
241
242    LOGICAL ::  switch_off_module = .FALSE.  !< local namelist parameter to switch off the module
243                                             !< although the respective module namelist appears in
244                                             !< the namelist file
245
246    NAMELIST /virtual_flight_parameters/  flight_angle,                                            &
247                                          flight_begin,                                            &
248                                          flight_end,                                              &
249                                          flight_level,                                            &
250                                          leg_mode,                                                &
251                                          max_elev_change,                                         &
252                                          rate_of_climb,                                           &
253                                          speed_agl,                                               &
254                                          switch_off_module,                                       &
255                                          x_end,                                                   &
256                                          x_start,                                                 &
257                                          y_end,                                                   &
258                                          y_start
259
260
261!
262!-- Move to the beginning of the namelist file and try to find and read the namelist.
263    REWIND( 11 )
264    READ( 11, virtual_flight_parameters, IOSTAT=io_status )
265
266!
267!-- Action depending on the READ status
268    IF ( io_status == 0 )  THEN
269!
270!--    virtual_flight_parameters namelist was found and read correctly. Set switch that virtual
271!--    flights are carried out.
272       IF ( .NOT. switch_off_module )  virtual_flight = .TRUE.
273
274    ELSEIF ( io_status > 0 )  THEN
275!
276!--    virtual_flight_parameters namelist was found, but contained errors. Print an error message
277!--    including the line that caused the problem.
278       BACKSPACE( 11 )
279       READ( 11 , '(A)' ) line
280       CALL parin_fail_message( 'virtual_flight_parameters', line )
281
282    ENDIF
283
284 END SUBROUTINE flight_parin
285
286
287!--------------------------------------------------------------------------------------------------!
288! Description:
289! ------------
290!> Inititalization of required arrays, number of legs and flags. Moreover, initialize flight speed
291!> and -direction, as well as start positions.
292!--------------------------------------------------------------------------------------------------!
293 SUBROUTINE flight_init
294
295    USE basic_constants_and_equations_mod,                                                         &
296        ONLY:  pi
297
298    USE control_parameters,                                                                        &
299        ONLY:  initializing_actions
300
301    USE indices,                                                                                   &
302        ONLY:  nxlg,                                                                               &
303               nxrg,                                                                               &
304               nysg,                                                                               &
305               nyng,                                                                               &
306               nzb,                                                                                &
307               nzt
308
309    IMPLICIT NONE
310
311    REAL(wp) ::  distance  !< distance between start and end position of a flight leg
312
313
314    IF ( debug_output )  CALL debug_message( 'flight_init', 'start' )
315!
316!-- Determine the number of flight legs
317    l = 1
318    DO  WHILE ( x_start(l) /= 999999999.0_wp  .AND.  l <= SIZE(x_start) )
319       l       = l + 1
320    ENDDO
321    num_leg = l-1
322!
323!-- Check for proper parameter settings
324    CALL flight_check_parameters
325!
326!-- Allocate and initialize logical array for flight pattern
327    ALLOCATE( cyclic_leg(1:num_leg) )
328!
329!-- Initialize flags for cyclic/return legs
330    DO  l = 1, num_leg
331       cyclic_leg(l) = MERGE( .TRUE., .FALSE., TRIM( leg_mode(l) ) == 'cyclic' )
332    ENDDO
333!
334!-- Allocate and initialize arraxs for flight position and speed. In case of restart runs these data
335!-- are read by the routine read_flight_restart_data instead.
336    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
337
338       ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
339!
340!--    Inititalize x-, y-, and z-positions with initial start position
341       x_pos(1:num_leg) = x_start(1:num_leg)
342       y_pos(1:num_leg) = y_start(1:num_leg)
343       z_pos(1:num_leg) = flight_level(1:num_leg)
344!
345!--    Allocate arrays for flight-speed components
346       ALLOCATE( u_agl(1:num_leg),                                                                 &
347                 v_agl(1:num_leg),                                                                 &
348                 w_agl(1:num_leg) )
349!
350!--    Inititalize u-, v- and w-component.
351       DO  l = 1, num_leg
352!
353!--       In case of return-legs, the flight direction, i.e. the horizontal flight-speed components,
354!--       are derived from the given start/end positions.
355          IF (  .NOT.  cyclic_leg(l) )  THEN
356             distance = SQRT( ( x_end(l) - x_start(l) )**2 + ( y_end(l) - y_start(l) )**2 )
357             u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
358             v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
359             w_agl(l) = rate_of_climb(l)
360!
361!--       In case of cyclic-legs, flight direction is directly derived from the given flight angle.
362          ELSE
363             u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
364             v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
365             w_agl(l) = rate_of_climb(l)
366          ENDIF
367
368       ENDDO
369
370    ENDIF
371!
372!-- Initialized data output
373    CALL flight_init_output
374!
375!-- Allocate array required for user-defined quantities if necessary.
376    IF ( num_var_fl_user  > 0 )  ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
377!
378!-- Allocate and initialize arrays containing the measured data
379    ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
380    ALLOCATE( sensor(1:num_var_fl,1:num_leg)   )
381    sensor_l = 0.0_wp
382    sensor   = 0.0_wp
383
384    IF ( debug_output )  CALL debug_message( 'flight_init', 'end' )
385
386 END SUBROUTINE flight_init
387
388
389!--------------------------------------------------------------------------------------------------!
390! Description:
391! ------------
392!> Initialization of output-variable names and units.
393!--------------------------------------------------------------------------------------------------!
394 SUBROUTINE flight_init_output
395
396    USE control_parameters,                                                                        &
397       ONLY:  cloud_droplets,                                                                      &
398              humidity,                                                                            &
399              neutral,                                                                             &
400              passive_scalar
401
402    USE bulk_cloud_model_mod,                                                                      &
403        ONLY:  bulk_cloud_model
404
405    USE netcdf_interface
406
407    IMPLICIT NONE
408
409    CHARACTER(LEN=10) ::  label_leg  !< dummy argument to convert integer to string
410
411    INTEGER(iwp) ::  i         !< loop variable
412    INTEGER(iwp) ::  id_pt     !< identifyer for labeling
413    INTEGER(iwp) ::  id_q      !< identifyer for labeling
414    INTEGER(iwp) ::  id_ql     !< identifyer for labeling
415    INTEGER(iwp) ::  id_s      !< identifyer for labeling
416    INTEGER(iwp) ::  id_u = 1  !< identifyer for labeling
417    INTEGER(iwp) ::  id_v = 2  !< identifyer for labeling
418    INTEGER(iwp) ::  id_w = 3  !< identifyer for labeling
419    INTEGER(iwp) ::  k         !< index variable
420
421    LOGICAL      ::  init = .TRUE.  !< flag to distiquish calls of user_init_flight
422!
423!-- Define output quanities, at least three variables are measured (u,v,w)
424    num_var_fl = 3
425    IF ( .NOT. neutral                     )  THEN
426       num_var_fl = num_var_fl + 1
427       id_pt      = num_var_fl
428    ENDIF
429    IF ( humidity                          )  THEN
430       num_var_fl = num_var_fl + 1
431       id_q       = num_var_fl
432    ENDIF
433    IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
434       num_var_fl = num_var_fl + 1
435       id_ql      = num_var_fl
436    ENDIF
437    IF ( passive_scalar                    )  THEN
438       num_var_fl = num_var_fl + 1
439       id_s       = num_var_fl
440    ENDIF
441!
442!-- Write output strings for dimensions x, y, z
443    DO  l=1, num_leg
444
445       IF ( l < 10                    )  WRITE( label_leg, '(I1)' )  l
446       IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)' )  l
447       IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)' )  l
448
449       dofl_dim_label_x(l)  = 'x_' // TRIM( label_leg )
450       dofl_dim_label_y(l)  = 'y_' // TRIM( label_leg )
451       dofl_dim_label_z(l)  = 'z_' // TRIM( label_leg )
452
453    ENDDO
454
455!
456!-- Call user routine to initialize further variables
457    CALL user_init_flight( init )
458!
459!-- Write output labels and units for the quanities
460    k = 1
461    DO  l=1, num_leg
462
463       IF ( l < 10                    )  WRITE( label_leg, '(I1)' )  l
464       IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)' )  l
465       IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)' )  l
466
467       label_leg = 'leg_' // TRIM(label_leg)
468       DO  i=1, num_var_fl
469
470          IF ( i == id_u      )  THEN
471             dofl_label(k) = TRIM( label_leg ) // '_u'
472             dofl_unit(k)  = 'm/s'
473             k             = k + 1
474          ELSEIF ( i == id_v  )  THEN
475             dofl_label(k) = TRIM( label_leg ) // '_v'
476             dofl_unit(k)  = 'm/s'
477             k             = k + 1
478          ELSEIF ( i == id_w  )  THEN
479             dofl_label(k) = TRIM( label_leg ) // '_w'
480             dofl_unit(k)  = 'm/s'
481             k             = k + 1
482          ELSEIF ( i == id_pt )  THEN
483             dofl_label(k) = TRIM( label_leg ) // '_theta'
484             dofl_unit(k)  = 'K'
485             k             = k + 1
486          ELSEIF ( i == id_q  )  THEN
487             dofl_label(k) = TRIM( label_leg ) // '_q'
488             dofl_unit(k)  = 'kg/kg'
489             k             = k + 1
490          ELSEIF ( i == id_ql )  THEN
491             dofl_label(k) = TRIM( label_leg ) // '_ql'
492             dofl_unit(k)  = 'kg/kg'
493             k             = k + 1
494          ELSEIF ( i == id_s  )  THEN
495             dofl_label(k) = TRIM( label_leg ) // '_s'
496             dofl_unit(k)  = 'kg/kg'
497             k             = k + 1
498          ENDIF
499       ENDDO
500
501       DO i=1, num_var_fl_user
502          CALL user_init_flight( init, k, i, label_leg )
503       ENDDO
504
505    ENDDO
506!
507!-- Finally, set the total number of flight-output quantities.
508    num_var_fl = num_var_fl + num_var_fl_user
509
510 END SUBROUTINE flight_init_output
511
512!--------------------------------------------------------------------------------------------------!
513! Description:
514! ------------
515!> This routine calculates the current flight positions and calls the respective interpolation
516!> routine to measure the data at the current flight position.
517!--------------------------------------------------------------------------------------------------!
518 SUBROUTINE flight_measurement
519
520    USE arrays_3d,                                                                                 &
521        ONLY:  ddzu,                                                                               &
522               ddzw,                                                                               &
523               pt,                                                                                 &
524               q,                                                                                  &
525               ql,                                                                                 &
526               s,                                                                                  &
527               u,                                                                                  &
528               v,                                                                                  &
529               w,                                                                                  &
530               zu,                                                                                 &
531               zw
532
533    USE control_parameters,                                                                        &
534        ONLY:  cloud_droplets,                                                                     &
535               dt_3d,                                                                              &
536               humidity,                                                                           &
537               neutral,                                                                            &
538               passive_scalar,                                                                     &
539               time_since_reference_point
540
541    USE cpulog,                                                                                    &
542        ONLY:  cpu_log,                                                                            &
543               log_point
544
545    USE grid_variables,                                                                            &
546        ONLY:  ddx,                                                                                &
547               ddy,                                                                                &
548               dx,                                                                                 &
549               dy
550
551    USE indices,                                                                                   &
552        ONLY:  nx,                                                                                 &
553               nxl,                                                                                &
554               nxr,                                                                                &
555               ny,                                                                                 &
556               nys,                                                                                &
557               nyn
558
559    USE bulk_cloud_model_mod,                                                                      &
560        ONLY:  bulk_cloud_model
561
562    USE pegrid
563
564    IMPLICIT NONE
565
566    INTEGER(iwp) ::  i  !< index of current grid box along x
567    INTEGER(iwp) ::  j  !< index of current grid box along y
568    INTEGER(iwp) ::  n  !< loop variable for number of user-defined output quantities
569
570    LOGICAL  ::  on_pe  !< flag to check if current flight position is on current PE
571
572    REAL(wp) ::  x  !< distance between left edge of current grid box and flight position
573    REAL(wp) ::  y  !< distance between south edge of current grid box and flight position
574
575    CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
576!
577!-- Perform flight measurement if start time is reached.
578    IF ( time_since_reference_point >= flight_begin  .AND.  time_since_reference_point <= flight_end )  THEN
579
580       sensor_l = 0.0_wp
581       sensor   = 0.0_wp
582!
583!--    Loop over all flight legs
584       DO  l = 1, num_leg
585!
586!--       Update location for each leg
587          x_pos(l) = x_pos(l) + u_agl(l) * dt_3d
588          y_pos(l) = y_pos(l) + v_agl(l) * dt_3d
589          z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
590!
591!--       Check if location must be modified for return legs.
592!--       Carry out horizontal reflection if required.
593          IF ( .NOT. cyclic_leg(l) )  THEN
594
595             IF ( x_start(l) <= x_end(l) )  THEN
596!
597!--             Outward flight, i.e. from start to end
598                IF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_end(l)      )  THEN
599                   x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
600                   u_agl(l) = - u_agl(l)
601!
602!--             Return flight, i.e. from end to start
603                ELSEIF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_start(l) )  THEN
604                   x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
605                   u_agl(l) = - u_agl(l)
606                ENDIF
607             ELSE
608!
609!--             Outward flight, i.e. from start to end
610                IF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_end(l)      )  THEN
611                   x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
612                   u_agl(l) = - u_agl(l)
613!
614!--             Return flight, i.e. from end to start
615                ELSEIF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_start(l) )  THEN
616                   x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
617                   u_agl(l) = - u_agl(l)
618                ENDIF
619             ENDIF
620
621             IF ( y_start(l) <= y_end(l) )  THEN
622!
623!--             Outward flight, i.e. from start to end
624                IF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_end(l)      )  THEN
625                   y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
626                   v_agl(l) = - v_agl(l)
627!
628!--             Return flight, i.e. from end to start
629                ELSEIF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_start(l) )  THEN
630                   y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
631                   v_agl(l) = - v_agl(l)
632                ENDIF
633             ELSE
634!
635!--             Outward flight, i.e. from start to end
636                IF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_end(l)      )  THEN
637                   y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
638                   v_agl(l) = - v_agl(l)
639!
640!--             Return flight, i.e. from end to start
641                ELSEIF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_start(l) )  THEN
642                   y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
643                   v_agl(l) = - v_agl(l)
644                ENDIF
645             ENDIF
646!
647!--       Check if flight position is outside the model domain and apply cyclic conditions if required
648          ELSEIF ( cyclic_leg(l) )  THEN
649!
650!--          Check if aircraft leaves the model domain at the right boundary
651             IF ( ( flight_angle(l) >= 0.0_wp     .AND.                                            &
652                    flight_angle(l) <= 90.0_wp )  .OR.                                             &
653                  ( flight_angle(l) >= 270.0_wp   .AND.                                            &
654                    flight_angle(l) <= 360.0_wp ) )  THEN
655                IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx )                                            &
656                     x_pos(l) = x_pos(l) - ( nx + 1 ) * dx
657!
658!--          Check if aircraft leaves the model domain at the left boundary
659             ELSEIF ( flight_angle(l) > 90.0_wp  .AND.  flight_angle(l) < 270.0_wp )  THEN
660                IF ( x_pos(l) < -0.5_wp * dx )                                                     &
661                     x_pos(l) = ( nx + 1 ) * dx + x_pos(l)
662             ENDIF
663!
664!--          Check if aircraft leaves the model domain at the north boundary
665             IF ( flight_angle(l) >= 0.0_wp  .AND.  flight_angle(l) <= 180.0_wp )  THEN
666                IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy )                                            &
667                     y_pos(l) = y_pos(l) - ( ny + 1 ) * dy
668!
669!--          Check if aircraft leaves the model domain at the south boundary
670             ELSEIF ( flight_angle(l) > 180.0_wp  .AND.  flight_angle(l) < 360.0_wp )  THEN
671                IF ( y_pos(l) < -0.5_wp * dy )                                                     &
672                     y_pos(l) = ( ny + 1 ) * dy + y_pos(l)
673             ENDIF
674
675          ENDIF
676!
677!--       Check if maximum elevation change is already reached. If required reflect vertically.
678          IF ( rate_of_climb(l) /= 0.0_wp )  THEN
679!
680!--          First check if aircraft is too high
681             IF (  w_agl(l) > 0.0_wp  .AND.  z_pos(l) - flight_level(l) > max_elev_change(l) )  THEN
682                z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) ) - z_pos(l)
683                w_agl(l) = - w_agl(l)
684!
685!--          Check if aircraft is too low
686             ELSEIF (  w_agl(l) < 0.0_wp  .AND.  z_pos(l) < flight_level(l) )  THEN
687                z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
688                w_agl(l) = - w_agl(l)
689             ENDIF
690
691          ENDIF
692!
693!--       Determine grid indices for flight position along x- and y-direction. Please note, there is
694!--       a special treatment for the index along z-direction, which is due to vertical grid stretching.
695          i = ( x_pos(l) + 0.5_wp * dx ) * ddx
696          j = ( y_pos(l) + 0.5_wp * dy ) * ddy
697!
698!--       Check if indices are on current PE
699          on_pe = ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )
700
701          IF ( on_pe )  THEN
702
703             var_index = 1
704!
705!--          Recalculate indices, as u is shifted by -0.5*dx.
706             i =   x_pos(l) * ddx
707             j = ( y_pos(l) + 0.5_wp * dy ) * ddy
708!
709!--          Calculate distance from left and south grid-box edges.
710             x  = x_pos(l) - ( 0.5_wp - i ) * dx
711             y  = y_pos(l) - j * dy
712!
713!--          Interpolate u-component onto current flight position.
714             CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
715             var_index = var_index + 1
716!
717!--          Recalculate indices, as v is shifted by -0.5*dy.
718             i = ( x_pos(l) + 0.5_wp * dx ) * ddx
719             j =   y_pos(l) * ddy
720
721             x  = x_pos(l) - i * dx
722             y  = y_pos(l) - ( 0.5_wp - j ) * dy
723             CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
724             var_index = var_index + 1
725!
726!--          Interpolate w and scalar quantities. Recalculate indices.
727             i  = ( x_pos(l) + 0.5_wp * dx ) * ddx
728             j  = ( y_pos(l) + 0.5_wp * dy ) * ddy
729             x  = x_pos(l) - i * dx
730             y  = y_pos(l) - j * dy
731!
732!--          Interpolate w-velocity component.
733             CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
734             var_index = var_index + 1
735!
736!--          Potential temerature
737             IF ( .NOT. neutral )  THEN
738                CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
739                var_index = var_index + 1
740             ENDIF
741!
742!--          Humidity
743             IF ( humidity )  THEN
744                CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
745                var_index = var_index + 1
746             ENDIF
747!
748!--          Liquid water content
749             IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
750                CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
751                var_index = var_index + 1
752             ENDIF
753!
754!--          Passive scalar
755             IF ( passive_scalar )  THEN
756                CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
757                var_index = var_index + 1
758             ENDIF
759!
760!--          Treat user-defined variables if required
761             DO  n = 1, num_var_fl_user
762                CALL user_flight( var_u, n )
763                CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
764                var_index = var_index + 1
765             ENDDO
766          ENDIF
767
768       ENDDO
769!
770!--    Write local data on global array.
771#if defined( __parallel )
772       CALL MPI_ALLREDUCE( sensor_l(1,1), sensor(1,1), num_var_fl * num_leg, MPI_REAL, MPI_SUM,    &
773                           comm2d, ierr )
774#else
775       sensor     = sensor_l
776#endif
777    ENDIF
778
779    CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )
780
781 END SUBROUTINE flight_measurement
782
783!--------------------------------------------------------------------------------------------------!
784! Description:
785! ------------
786!> This subroutine bi-linearly interpolates the respective data onto the current flight position.
787!--------------------------------------------------------------------------------------------------!
788 SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )
789
790    USE control_parameters,                                                                        &
791        ONLY:  dz,                                                                                 &
792               dz_stretch_level_start
793
794    USE grid_variables,                                                                            &
795       ONLY:  dx,                                                                                  &
796              dy
797
798    USE indices,                                                                                   &
799        ONLY:  nzb,                                                                                &
800               nzt,                                                                                &
801               nxlg,                                                                               &
802               nxrg,                                                                               &
803               nysg,                                                                               &
804               nyng
805
806    IMPLICIT NONE
807
808    INTEGER(iwp) ::  i        !< index along x
809    INTEGER(iwp) ::  j        !< index along y
810    INTEGER(iwp) ::  k        !< index along z
811    INTEGER(iwp) ::  k1       !< dummy variable
812    INTEGER(iwp) ::  var_ind  !< index variable for output quantity
813
814    REAL(wp) ::  aa         !< dummy argument for horizontal interpolation
815    REAL(wp) ::  bb         !< dummy argument for horizontal interpolation
816    REAL(wp) ::  cc         !< dummy argument for horizontal interpolation
817    REAL(wp) ::  dd         !< dummy argument for horizontal interpolation
818    REAL(wp) ::  gg         !< dummy argument for horizontal interpolation
819    REAL(wp) ::  fac        !< flag to indentify if quantity is on zu or zw level
820    REAL(wp) ::  var_int    !< horizontal interpolated variable at current position
821    REAL(wp) ::  var_int_l  !< horizontal interpolated variable at k-level
822    REAL(wp) ::  var_int_u  !< horizontal interpolated variable at (k+1)-level
823    REAL(wp) ::  x          !< distance between left edge of current grid box and flight position
824    REAL(wp) ::  y          !< distance between south edge of current grid box and flight position
825
826    REAL(wp), DIMENSION(1:nzt+1)   ::  ddz_uw  !< inverse vertical grid spacing
827    REAL(wp), DIMENSION(nzb:nzt+1) ::  z_uw    !< height level
828
829    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< treated quantity
830!
831!-- Calculate interpolation coefficients
832    aa = x**2          + y**2
833    bb = ( dx - x )**2 + y**2
834    cc = x**2          + ( dy - y )**2
835    dd = ( dx - x )**2 + ( dy - y )**2
836    gg = aa + bb + cc + dd
837!
838!-- Obtain vertical index. Special treatment for grid index along z-direction if flight position is
839!-- above the vertical grid-stretching level. fac=1 if variable is on scalar grid level, fac=0 for
840!-- w-component.
841    IF ( z_pos(l) < dz_stretch_level_start(1) )  THEN
842       k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
843    ELSE
844!
845!--    Search for k-index
846       DO  k1 = nzb, nzt
847          IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) )  THEN
848             k = k1
849             EXIT
850          ENDIF
851       ENDDO
852    ENDIF
853!
854!-- (x,y)-interpolation at k-level
855    var_int_l = ( ( gg - aa ) * var(k,j,i)       +                                                 &
856                  ( gg - bb ) * var(k,j,i+1)     +                                                 &
857                  ( gg - cc ) * var(k,j+1,i)     +                                                 &
858                  ( gg - dd ) * var(k,j+1,i+1)                                                     &
859                ) / ( 3.0_wp * gg )
860!
861!-- (x,y)-interpolation on (k+1)-level
862    var_int_u = ( ( gg - aa ) * var(k+1,j,i)     +                                                 &
863                  ( gg - bb ) * var(k+1,j,i+1)   +                                                 &
864                  ( gg - cc ) * var(k+1,j+1,i)   +                                                 &
865                  ( gg - dd ) * var(k+1,j+1,i+1)                                                   &
866                ) / ( 3.0_wp * gg )
867!
868!-- z-interpolation onto current flight postion
869    var_int = var_int_l + ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1) * (var_int_u - var_int_l )
870!
871!-- Store on local data array
872    sensor_l(var_ind,l) = var_int
873
874 END SUBROUTINE interpolate_xyz
875
876!--------------------------------------------------------------------------------------------------!
877! Description:
878! ------------
879!> Perform parameter checks.
880!--------------------------------------------------------------------------------------------------!
881 SUBROUTINE flight_check_parameters
882
883    USE arrays_3d,                                                                                 &
884        ONLY:  zu
885
886    USE control_parameters,                                                                        &
887        ONLY:  bc_lr_cyc,                                                                          &
888               bc_ns_cyc,                                                                          &
889               message_string
890
891    USE grid_variables,                                                                            &
892        ONLY:  dx,                                                                                 &
893               dy
894
895    USE indices,                                                                                   &
896        ONLY:  nx,                                                                                 &
897               ny,                                                                                 &
898               nz
899
900    USE netcdf_interface,                                                                          &
901        ONLY:  netcdf_data_format
902
903    IMPLICIT NONE
904
905!
906!-- Check if start positions are properly set.
907    DO  l = 1, num_leg
908       IF ( x_start(l) < 0.0_wp  .OR.  x_start(l) > ( nx + 1 ) * dx )  THEN
909          message_string = 'Start x position is outside the model domain'
910          CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
911       ENDIF
912       IF ( y_start(l) < 0.0_wp  .OR.  y_start(l) > ( ny + 1 ) * dy )  THEN
913          message_string = 'Start y position is outside the model domain'
914          CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
915       ENDIF
916
917    ENDDO
918!
919!-- Check for leg mode
920    DO  l = 1, num_leg
921!
922!--    Check if leg mode matches the overall lateral model boundary conditions.
923       IF ( TRIM( leg_mode(l) ) == 'cyclic' )  THEN
924          IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
925             message_string = 'Cyclic flight leg does not match lateral boundary condition'
926             CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
927          ENDIF
928!
929!--    Check if end-positions are inside the model domain in case of return-legs.
930       ELSEIF ( TRIM( leg_mode(l) ) == 'return' )  THEN
931          IF ( x_end(l) > ( nx + 1 ) * dx  .OR.  y_end(l) > ( ny + 1 ) * dx )  THEN
932             message_string = 'Flight leg or parts of it are outside the model domain'
933             CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
934          ENDIF
935       ELSE
936          message_string = 'Unknown flight mode'
937          CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
938       ENDIF
939
940    ENDDO
941!
942!-- Check if given flight object remains inside model domain if a rate of climb / descent is
943!-- prescribed.
944    DO  l = 1, num_leg
945       IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR.                                     &
946            flight_level(l) + max_elev_change(l) <= 0.0_wp )  THEN
947          message_string = 'Flight level is outside the model domain '
948          CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
949       ENDIF
950    ENDDO
951!
952!-- Check for appropriate NetCDF format. Definition of more than one unlimited dimension is
953!-- unfortunately only possible with NetCDF4/HDF5.
954    IF (  netcdf_data_format <= 2 )  THEN
955       message_string = 'netcdf_data_format must be > 2'
956       CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
957    ENDIF
958
959
960 END SUBROUTINE flight_check_parameters
961
962
963!--------------------------------------------------------------------------------------------------!
964! Description:
965! ------------
966!> Read module-specific global restart data (Fortran binary format).
967!--------------------------------------------------------------------------------------------------!
968 SUBROUTINE flight_rrd_global_ftn( found )
969
970
971    USE control_parameters,                                                                        &
972        ONLY: length,                                                                              &
973              restart_string
974
975
976    IMPLICIT NONE
977
978    LOGICAL, INTENT(OUT)  ::  found  !< flag indicating if a variable string is found
979
980
981    found = .TRUE.
982
983
984    SELECT CASE ( restart_string(1:length) )
985
986       CASE ( 'u_agl' )
987          IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
988          READ ( 13 )  u_agl
989       CASE ( 'v_agl' )
990          IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
991          READ ( 13 )  v_agl
992       CASE ( 'w_agl' )
993          IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
994          READ ( 13 )  w_agl
995       CASE ( 'x_pos' )
996          IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
997          READ ( 13 )  x_pos
998       CASE ( 'y_pos' )
999          IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
1000          READ ( 13 )  y_pos
1001       CASE ( 'z_pos' )
1002          IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
1003          READ ( 13 )  z_pos
1004
1005       CASE DEFAULT
1006
1007          found = .FALSE.
1008
1009    END SELECT
1010
1011
1012 END SUBROUTINE flight_rrd_global_ftn
1013
1014
1015 
1016!------------------------------------------------------------------------------!
1017! Description:
1018! ------------
1019!> Read module-specific global restart data (MPI-IO).
1020!------------------------------------------------------------------------------!
1021 SUBROUTINE flight_rrd_global_mpi
1022
1023
1024    IMPLICIT NONE
1025
1026    LOGICAL  ::  array_found  !< flag indicating if respective array is found in restart file
1027
1028
1029    CALL rd_mpi_io_check_array( 'u_agl', found = array_found )
1030    IF ( array_found)  THEN
1031       IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
1032       CALL rrd_mpi_io_global_array( 'u_agl', u_agl )
1033    ENDIF
1034    CALL rd_mpi_io_check_array( 'v_agl', found = array_found )
1035    IF ( array_found)  THEN
1036       IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
1037       CALL rrd_mpi_io_global_array( 'v_agl', v_agl )
1038    ENDIF
1039    CALL rd_mpi_io_check_array( 'w_agl', found = array_found )
1040    IF ( array_found)  THEN
1041       IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
1042       CALL rrd_mpi_io_global_array( 'w_agl', w_agl )
1043    ENDIF
1044    CALL rd_mpi_io_check_array( 'x_pos', found = array_found )
1045    IF ( array_found)  THEN
1046       IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
1047       CALL rrd_mpi_io_global_array( 'x_pos', x_pos )
1048    ENDIF
1049    CALL rd_mpi_io_check_array( 'y_pos', found = array_found )
1050    IF ( array_found)  THEN
1051       IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
1052       CALL rrd_mpi_io_global_array( 'y_pos', y_pos )
1053    ENDIF
1054    CALL rd_mpi_io_check_array( 'z_pos', found = array_found )
1055    IF ( array_found)  THEN
1056       IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
1057       CALL rrd_mpi_io_global_array( 'z_pos', z_pos )
1058    ENDIF
1059
1060 END SUBROUTINE flight_rrd_global_mpi
1061
1062   
1063    !--------------------------------------------------------------------------------------------------!
1064! Description:
1065! ------------
1066!> This routine writes the respective restart data.
1067!--------------------------------------------------------------------------------------------------!
1068 SUBROUTINE flight_wrd_global
1069
1070
1071    IMPLICIT NONE
1072
1073
1074    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
1075
1076       CALL wrd_write_string( 'u_agl' )
1077       WRITE ( 14 )  u_agl
1078
1079       CALL wrd_write_string( 'v_agl' )
1080       WRITE ( 14 )  v_agl
1081
1082       CALL wrd_write_string( 'w_agl' )
1083       WRITE ( 14 )  w_agl
1084
1085       CALL wrd_write_string( 'x_pos' )
1086       WRITE ( 14 )  x_pos
1087
1088       CALL wrd_write_string( 'y_pos' )
1089       WRITE ( 14 )  y_pos
1090
1091       CALL wrd_write_string( 'z_pos' )
1092       WRITE ( 14 )  z_pos
1093
1094    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
1095
1096       CALL wrd_mpi_io_global_array( 'u_agl', u_agl )
1097       CALL wrd_mpi_io_global_array( 'v_agl', v_agl )
1098       CALL wrd_mpi_io_global_array( 'w_agl', w_agl )
1099       CALL wrd_mpi_io_global_array( 'x_pos', x_pos )
1100       CALL wrd_mpi_io_global_array( 'y_pos', y_pos )
1101       CALL wrd_mpi_io_global_array( 'z_pos', z_pos )
1102
1103    ENDIF
1104
1105 END SUBROUTINE flight_wrd_global
1106
1107
1108 END MODULE flight_mod
Note: See TracBrowser for help on using the repository browser.