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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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