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

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