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

Last change on this file since 3053 was 3049, checked in by Giersch, 7 years ago

Revision history corrected

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