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

Last change on this file since 3173 was 3065, checked in by Giersch, 6 years ago

New vertical stretching procedure has been introduced

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