source: palm/trunk/SOURCE/virtual_flight_mod.f90 @ 4535

Last change on this file since 4535 was 4535, checked in by raasch, 12 months ago

bugfix for restart data format query

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