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

Last change on this file since 2807 was 2776, checked in by Giersch, 6 years ago

Skipping of module related restart data changed + adapting synthetic turbulence generator to current restart procedure

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