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

Last change on this file since 3899 was 3885, checked in by kanani, 5 years ago

restructure/add location/debug messages

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