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

Last change on this file since 2225 was 2101, checked in by suehring, 7 years ago

last commit documented

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