3.5.4
User-defined output quantities
A very typical request of users is the
calculation and
output of
quantities which are not part of PALM's standard output. The basic user
interface includes a number of subroutines which allow the calculation
of user-defined quantities and output of these quantities as time series, (horizontally averaged) vertical profile, 2d cross
section or 3d volume data. The respective subroutines
contain sample code lines (written as comment lines) for defining, calculating and
output of such quantities.
Output times, averaging intervals, etc. are steered by the same variables as used for the standard PALM output quantities, e.g. dt_data_output.
The
rest of this chapter explains step-by-step how to modify/extend the
default file user_interface.f90 in order to generate the respective output.
Output of vertical profiles
This example shows the output of the
quantity "turbulent resolved-scale horizontal momentum flux" (u*v*). If more than one user-defined
quantity shall be output, the following steps have to be carried out in the
same way for each of the quantities.
- The
quantity has to be given a unique string identifier, e.g. 'u*v*'.
This identifier must be different from the identifiers used for the
PALM standard output (see list in description of parameter data_output_pr).
To switch on output of this quantity, the user has to assign the string
identifier to the parameter data_output_pr_user,
eg.:
data_output_pr_user = 'u*v*',
- For
the
quantity, an identification number, a physical unit, and the vertical
grid on which it is defined (u- or w-grid), has to be assigned (subroutine user_check_data_output_pr):
CASE (
'u*v*' )
index = pr_palm + 1
! identification number
dopr_index(var_count) = index
dopr_unit(var_count) = 'm2/s2'
! physical unit
hom(:,2,index,:) = SPREAD( zu, 2, statistic_regions+1 ) ! vertical grid
Here only the those parts in red color have to be given by the user appropriately.
The identification number (index) must be within the range [ pr_palm+1 , pr_palm+max_pr_user ], where max_pr_user is the number of user-defined profiles as given by parameter data_output_pr_user in the respective PALM run. The physical unit has to be given with respect to the NetCDF conventions. If no unit is given,
PALM will abort. The vertical grid has to be either zu (u-grid) or zw (w-grid).
- The quantity has to be calculated for all gridpoints (subroutine user_statistics):
!$OMP DO
DO i = nxl, nxr
DO j = nys, nyn
DO k = nzb_s_outer(j,i)+1, nzt
sums_l(k,pr_palm+1,tn) = sums_l(k,pr_palm+1,tn) + &
( 0.5 * ( u(k,j,i) + u(k,j,i+1) ) - hom(k,1,1,sr) ) * &
( 0.5 * ( v(k,j,i) + v(k,j+1,i) ) - hom(k,1,2,sr) ) * &
* rmask(j,i,sr)
ENDDO
ENDDO
ENDDO
Once again, only those parts in red have to be adjusted by the user.
The
turbulent resolved-scale momentum flux u*v* is defined as the product
of the deviations of the horizontal velocities from their respective
horizontally averaged mean values. These mean values are stored in
array hom(..,1,1,sr) and hom(..,1,2,sr) for the u- and v-component, respectively. Since due to the staggered grid, u and v
are not defined at the same gridpoints, they have to be interpolated
appropriately (here to the center of the gridbox). The result of the
calculation is stored in array sums_l.
The second index of this array is the identification number of the
profile which must match the one given in the previous step 2.
Output of timeseries
still to be added
Output of 2d cross sections or 3d volume data
This example shows the output of the
quantity "square of the u-component" (Note: this quantity
could of course easily be calculated from the u-component by
postprocessing the PALM output so that calculation within PALM is not
necessarily required). If more than one user-defined
quantity shall be output, the following steps have to be carried out in the
same way for each of the quantities.
- The
quantity has to be given a unique string identifier, e.g. 'u2'.
This identifier must be different from the identifiers used for the
PALM standard output (see list in description of parameter data_output).
To switch on output of this quantity, the user has to assign the string
identifier to the parameter data_output_user,
eg.:
data_output_user = 'u2', 'u2_xy_av'
The
pure string 'u2'
switches on the output of instantaneous 3d volume data. Output of cross
section data and time averaged data is switched on by additionally
appending the strings '_xy',
'_xz', '_yz', and/or '_av' (for a
detailed explanation see parameter data_output).
- In
order to store the quantities' grid point data within PALM, a 3d data
array has to be declared in module user:
REAL,
DIMENSION(:,:,:), ALLOCATABLE :: u2, u2_av
The
second array u2_av
is needed in case that output of time averaged data is requested. It is
used to store the sum of the data of the respective time levels over
which the average has to be carried out.
- The
data array has to be allocated in subroutine user_init:
ALLOCATE( u2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
In
case that output of time averaged data is requested, the array
containing the sum has possibly to be read from the restart file (local
filename BININ)
by executing the following code in user_init:
IF ( initializing_actions == 'read_restart_data' )
THEN
READ ( 13 ) field_chr
DO WHILE ( TRIM( field_chr ) /= '*** end
user ***' )
SELECT CASE ( TRIM( field_chr ) )
CASE ( 'u2_av' )
ALLOCATE(
u2_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
READ ( 13 ) u2_av
CASE DEFAULT
PRINT*,
'+++ user_init: unknown
variable named "', &
TRIM(
field_chr ), '" found in'
PRINT*,
'
data from prior run on PE ', myid
CALL local_stop
END SELECT
ENDDO
ENDIF
- The
quantity has to be given a unit (subroutine user_check_data_output):
CASE (
'u2' )
unit = 'm2/s2'
Otherwise,
PALM will abort.
- The
vertical grid on which the quantity is defined (given by the levels
'zu' or 'zw', on which the u- or w-component of the velocity are
defined) has to be specified for the NetCDF output files in subroutine user_define_netcdf_grid:
CASE ( 'u2', 'u2_xy', 'u2_xz', 'u2_yz' )
grid = 'zu'
As
the example shows, this grid has to be defined for the 3d volume data
as well as for all of the three cross sections.
- After
each timestep, the quantity has to be calculated at all gridpoints and
to be stored. This has to be done in subroutine user_actions
at location 'after_integration':
CASE
( 'after_integration' )
!
!--
Enter actions to be done after every time integration (before
!--
data output)
!--
Sample for user-defined output:
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nzt+1
u2(k,j,i) = u(k,j,i)**2
ENDDO
ENDDO
ENDDO
- In
case that output of time-averaged data is requested, the sum- and
average-operations as well as the allocation of the sum-array have to
be carried out in subroutine user_3d_data_averaging:
IF ( mode == 'allocate' ) THEN
...
CASE ( 'u2' )
IF ( .NOT. ALLOCATED( u2_av ) ) THEN
ALLOCATE( u2_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
ENDIF
u2_av = 0.0
...
ELSEIF ( mode == 'sum' ) THEN
...
CASE
( 'u2' )
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nzt+1
u2_av(k,j,i) = u2_av(k,j,i) + u2(k,j,i)
ENDDO
ENDDO
ENDDO
...
ELSEIF ( mode == 'average' ) THEN
...
CASE
( 'u2' )
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nzt+1
u2_av(k,j,i) = u2_av(k,j,i) / REAL( average_count_3d )
ENDDO
ENDDO
ENDDO
- For
output of 2d cross sections, the gridpoint data of the quantity has to
be resorted to array local_pf
in subroutine user_data_output_2d.
Also the vertical grid, on which the quantity is defined, has to be set
again:
CASE
( 'u2_xy', 'u2_xz', 'u2_yz' )
IF ( av == 0 ) THEN
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nzt+1
local_pf(i,j,k) = u2(k,j,i)
ENDDO
ENDDO
ENDDO
ELSE
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nzt+1
local_pf(i,j,k) = u2_av(k,j,i)
ENDDO
ENDDO
ENDDO
ENDIF
grid = 'zu'
The ELSE case is
only needed in case that output of time-averaged data is requested.
- For
output of 3d volume data, the gridpoint data of the quantity has to be
resorted to array local_pf
in subroutine user_data_output_3d.
Also the vertical grid, on which the quantity is defined, has to be set
again:
CASE
( 'u2' )
IF ( av == 0 ) THEN
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nz_do
local_pf(i,j,k) = u2(k,j,i)
ENDDO
ENDDO
ENDDO
ELSE
DO i = nxl-1, nxr+1
DO j = nys-1, nyn+1
DO k = nzb, nz_do
local_pf(i,j,k) = u2_av(k,j,i)
ENDDO
ENDDO
ENDDO
ENDIF
grid = 'zu'
The ELSE case is
only needed in case that output of time-averaged data is requested.
- In
case of job chains, the sum array has to be written to the (binary)
restart file (local filename BINOUT)
in subroutine user_last_actions:
IF ( ALLOCATED( u2_av ) ) THEN
WRITE ( 14 )
'u2_av
'; WRITE ( 14 ) u2_av
ENDIF
Otherwise, the calculated
time-average may be wrong.
Last
change: $Id: chapter_3.5.4.html 89 2007-05-25 12:08:31Z heinze $