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.

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 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, these steps have to be carried out in the same way for each of the quantities.

  1. 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).

  2. 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.


  3. 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


  4. The quantity has to be given a unit (subroutine user_check_data_output):

        CASE ( 'u2' )
           unit = 'm2/s2'
     
    Otherwise, PALM will abort.


  5. 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.


  6. 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

     
  7. 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

     
  8. 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.


  9. 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.


  10. 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 87 2007-05-22 15:46:47Z raasch $