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

Last change on this file since 4004 was 4004, checked in by suehring, 2 years ago

chemistry: perform basic checks only when anthropenic emissions are switched on; virtual flights: allow arbitrary start/end positions also in return mode; bugfix in 2d data output

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