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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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