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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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