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

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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