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

Last change on this file since 1960 was 1960, checked in by suehring, 8 years ago

Separate balance equations for humidity and passive_scalar

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