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

Last change on this file since 3872 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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