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

Last change on this file since 1958 was 1958, checked in by suehring, 8 years ago

last commit documented

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