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

Last change on this file since 2932 was 2932, checked in by maronga, 3 years ago

renamed all Fortran NAMELISTS

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