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

Last change on this file since 2894 was 2894, checked in by Giersch, 3 years ago

Reading/Writing? data in case of restart runs revised

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