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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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