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

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

Forced header and separation lines into 80 columns

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