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

Last change on this file since 2563 was 2563, checked in by Giersch, 4 years ago

Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data

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