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

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

last commit documented

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