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

Last change on this file since 3540 was 3421, checked in by gronemeier, 6 years ago

new surface-data output; renamed output variables (pt to theta, rho_air to rho, rho_ocean to rho_sea_water)

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