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

Last change on this file since 4598 was 4541, checked in by suehring, 4 years ago

virtual flights: Bugfix, use time_since_reference_point instead of simulated_time

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