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

Last change on this file since 4495 was 4495, checked in by raasch, 4 years ago

restart data handling with MPI-IO added, first part

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