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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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