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

Last change on this file since 3253 was 3248, checked in by sward, 6 years ago

Minor format changes

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