Ignore:
Timestamp:
Jun 18, 2019 10:32:41 AM (5 years ago)
Author:
suehring
Message:

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90

    r3998 r4039  
    2525! -----------------
    2626! $Id$
     27! - Add output of uu, vv, ww to enable variance calculation according temporal
     28!   EC method
     29! - Allocate arrays only when they are required
     30! - Formatting adjustment
     31! - Rename subroutines
     32! - Further modularization
     33!
     34! 3998 2019-05-23 13:38:11Z suehring
    2735! Bugfix in gathering all output strings
    2836!
     
    4351 MODULE diagnostic_output_quantities_mod
    4452 
    45 
    46 !     USE arrays_3d,                                                             &
    47 !         ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
    48 !                precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
    49 !                u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
     53    USE arrays_3d,                                                             &
     54        ONLY:  ddzu, u, v, w
    5055!
    5156!     USE averaging
     
    6267!     USE cpulog,                                                                &
    6368!         ONLY:  cpu_log, log_point
    64 !
    65 !     USE indices,                                                               &
    66 !         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
    67 !                nzb, nzt, wall_flags_0
     69
     70   USE grid_variables,                                                         &
     71        ONLY:  ddx, ddy
     72!
     73    USE indices,                                                               &
     74        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
    6875!
    6976    USE kinds
     
    103110
    104111    INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
    105 
    106     LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE.
    107     LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.
     112 
     113    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
     114    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
    108115
    109116    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
    110117    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
     118   
     119    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu     !< uu
     120    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu
     121   
     122    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv     !< vv
     123    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv
     124   
     125    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
     126    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
    111127
    112128
     
    117133!
    118134!-- Public variables
    119     PUBLIC do_all,                                                                                 &
    120            initialized_diagnostic_output_quantities,                                               &
    121            prepared_diagnostic_output_quantities,                                                  &
    122            ti, ti_av,                                                                              &
    123            timestep_number_at_prev_calc
    124 !
    125 !-- Public routines
    126     PUBLIC diagnostic_output_quantities_calculate,                                                 &
    127            diagnostic_output_quantities_init,                                                      &
    128            diagnostic_output_quantities_prepare
    129 
    130 
    131     INTERFACE diagnostic_output_quantities_init
    132        MODULE PROCEDURE diagnostic_output_quantities_init
    133     END INTERFACE diagnostic_output_quantities_init
    134 
    135     INTERFACE diagnostic_output_quantities_calculate
    136        MODULE PROCEDURE diagnostic_output_quantities_calculate
    137     END INTERFACE diagnostic_output_quantities_calculate
    138 
    139     INTERFACE diagnostic_output_quantities_prepare
    140        MODULE PROCEDURE diagnostic_output_quantities_prepare
    141     END INTERFACE diagnostic_output_quantities_prepare
     135    PUBLIC do_all,                                                             &
     136           initialized_diagnostic_output_quantities,                           &
     137           prepared_diagnostic_output_quantities,                              &
     138           timestep_number_at_prev_calc,                                       &
     139           ti_av,                                                              &
     140           uu_av,                                                              &
     141           vv_av,                                                              &
     142           ww_av
     143!                                                                             
     144!-- Public routines                                                           
     145    PUBLIC doq_3d_data_averaging,                                              &
     146           doq_calculate,                                                      &
     147           doq_check_data_output,                                              &
     148           doq_define_netcdf_grid,                                             &
     149           doq_output_2d,                                                      &
     150           doq_output_3d,                                                      &
     151           doq_output_mask,                                                    &
     152           doq_init,                                                           &
     153           doq_prepare,                                                        &
     154           doq_wrd_local
     155!          doq_rrd_local,                                                      &
     156
     157
     158    INTERFACE doq_3d_data_averaging
     159       MODULE PROCEDURE doq_3d_data_averaging
     160    END INTERFACE doq_3d_data_averaging       
     161
     162    INTERFACE doq_calculate
     163       MODULE PROCEDURE doq_calculate
     164    END INTERFACE doq_calculate
     165
     166    INTERFACE doq_check_data_output
     167       MODULE PROCEDURE doq_check_data_output
     168    END INTERFACE doq_check_data_output
     169   
     170    INTERFACE doq_define_netcdf_grid
     171       MODULE PROCEDURE doq_define_netcdf_grid
     172    END INTERFACE doq_define_netcdf_grid
     173   
     174    INTERFACE doq_output_2d
     175       MODULE PROCEDURE doq_output_2d
     176    END INTERFACE doq_output_2d
     177   
     178    INTERFACE doq_output_3d
     179       MODULE PROCEDURE doq_output_3d
     180    END INTERFACE doq_output_3d
     181   
     182    INTERFACE doq_output_mask
     183       MODULE PROCEDURE doq_output_mask
     184    END INTERFACE doq_output_mask
     185     
     186    INTERFACE doq_init
     187       MODULE PROCEDURE doq_init
     188    END INTERFACE doq_init
     189
     190    INTERFACE doq_prepare
     191       MODULE PROCEDURE doq_prepare
     192    END INTERFACE doq_prepare
     193   
     194!     INTERFACE doq_rrd_local
     195!        MODULE PROCEDURE doq_rrd_local
     196!     END INTERFACE doq_rrd_local
     197   
     198    INTERFACE doq_wrd_local
     199       MODULE PROCEDURE doq_wrd_local
     200    END INTERFACE doq_wrd_local
    142201
    143202
    144203 CONTAINS
    145 
     204 
    146205!------------------------------------------------------------------------------!
    147206! Description:
    148207! ------------
    149 !> ...
    150 !------------------------------------------------------------------------------!
    151  SUBROUTINE diagnostic_output_quantities_init
    152 
    153 
    154     USE indices,                                                                                   &
    155         ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     208!> Sum up and time-average diagnostic output quantities as well as allocate
     209!> the array necessary for storing the average.
     210!------------------------------------------------------------------------------!
     211 SUBROUTINE doq_3d_data_averaging( mode, variable )
     212
     213    USE control_parameters,                                                    &
     214        ONLY:  average_count_3d
     215
     216    CHARACTER (LEN=*) ::  mode    !<
     217    CHARACTER (LEN=*) :: variable !<
     218
     219    INTEGER(iwp) ::  i !<
     220    INTEGER(iwp) ::  j !<
     221    INTEGER(iwp) ::  k !<
     222
     223    IF ( mode == 'allocate' )  THEN
     224
     225       SELECT CASE ( TRIM( variable ) )
     226
     227          CASE ( 'ti' )
     228             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
     229                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     230             ENDIF
     231             ti_av = 0.0_wp
     232       
     233          CASE ( 'uu' )
     234             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
     235                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     236             ENDIF
     237             uu_av = 0.0_wp
     238               
     239          CASE ( 'vv' )
     240             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
     241                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     242             ENDIF
     243             vv_av = 0.0_wp
     244               
     245          CASE ( 'ww' )
     246             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
     247                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     248             ENDIF
     249             ww_av = 0.0_wp
     250
     251          CASE DEFAULT
     252             CONTINUE
     253
     254       END SELECT
     255
     256    ELSEIF ( mode == 'sum' )  THEN
     257
     258       SELECT CASE ( TRIM( variable ) )
     259 
     260          CASE ( 'ti' )
     261             IF ( ALLOCATED( ti_av ) ) THEN
     262                DO  i = nxl, nxr
     263                   DO  j = nys, nyn
     264                      DO  k = nzb, nzt+1
     265                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
     266                      ENDDO
     267                   ENDDO
     268                ENDDO
     269             ENDIF
     270             
     271          CASE ( 'uu' )
     272             IF ( ALLOCATED( uu_av ) ) THEN
     273                DO  i = nxl, nxr
     274                   DO  j = nys, nyn
     275                      DO  k = nzb, nzt+1
     276                         uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
     277                      ENDDO
     278                   ENDDO
     279                ENDDO
     280             ENDIF
     281             
     282          CASE ( 'vv' )
     283             IF ( ALLOCATED( vv_av ) ) THEN
     284                DO  i = nxl, nxr
     285                   DO  j = nys, nyn
     286                      DO  k = nzb, nzt+1
     287                         vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
     288                      ENDDO
     289                   ENDDO
     290                ENDDO
     291             ENDIF
     292             
     293          CASE ( 'ww' )
     294             IF ( ALLOCATED( ww_av ) ) THEN
     295                DO  i = nxl, nxr
     296                   DO  j = nys, nyn
     297                      DO  k = nzb, nzt+1
     298                         ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
     299                      ENDDO
     300                   ENDDO
     301                ENDDO
     302             ENDIF
     303
     304          CASE DEFAULT
     305             CONTINUE
     306
     307       END SELECT
     308
     309    ELSEIF ( mode == 'average' )  THEN
     310
     311       SELECT CASE ( TRIM( variable ) )
     312
     313          CASE ( 'ti' )
     314             IF ( ALLOCATED( ti_av ) ) THEN
     315                DO  i = nxl, nxr
     316                   DO  j = nys, nyn
     317                      DO  k = nzb, nzt+1
     318                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     319                      ENDDO
     320                   ENDDO
     321                ENDDO
     322             ENDIF
     323       
     324          CASE ( 'uu' )
     325             IF ( ALLOCATED( uu_av ) ) THEN
     326                DO  i = nxl, nxr
     327                   DO  j = nys, nyn
     328                      DO  k = nzb, nzt+1
     329                         uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     330                      ENDDO
     331                   ENDDO
     332                ENDDO
     333             ENDIF
     334             
     335          CASE ( 'vv' )
     336             IF ( ALLOCATED( vv_av ) ) THEN
     337                DO  i = nxl, nxr
     338                   DO  j = nys, nyn
     339                      DO  k = nzb, nzt+1
     340                         vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     341                      ENDDO
     342                   ENDDO
     343                ENDDO
     344             ENDIF
     345             
     346          CASE ( 'ww' )
     347             IF ( ALLOCATED( ww_av ) ) THEN
     348                DO  i = nxl, nxr
     349                   DO  j = nys, nyn
     350                      DO  k = nzb, nzt+1
     351                         ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     352                      ENDDO
     353                   ENDDO
     354                ENDDO
     355             ENDIF
     356
     357       END SELECT
     358
     359    ENDIF
     360
     361
     362 END SUBROUTINE doq_3d_data_averaging
     363 
     364!------------------------------------------------------------------------------!
     365! Description:
     366! ------------
     367!> Check data output for diagnostic output
     368!------------------------------------------------------------------------------!
     369 SUBROUTINE doq_check_data_output( var, unit )
    156370
    157371    IMPLICIT NONE
     372
     373    CHARACTER (LEN=*) ::  unit  !<
     374    CHARACTER (LEN=*) ::  var   !<
     375
     376    SELECT CASE ( TRIM( var ) )
     377
     378       CASE ( 'ti' )
     379          unit = '1/s'
     380             
     381       CASE ( 'uu' )
     382          unit = 'm2/s2'
     383             
     384       CASE ( 'vv' )
     385          unit = 'm2/s2'
     386             
     387       CASE ( 'ww' )
     388          unit = 'm2/s2'
     389             
     390       CASE DEFAULT
     391          unit = 'illegal'
     392
     393    END SELECT
     394
     395
     396 END SUBROUTINE doq_check_data_output
     397 
     398!------------------------------------------------------------------------------!
     399!
     400! Description:
     401! ------------
     402!> Subroutine defining appropriate grid for netcdf variables.
     403!------------------------------------------------------------------------------!
     404 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
     405   
     406    IMPLICIT NONE
     407
     408    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
     409    LOGICAL, INTENT(OUT)           ::  found       !<
     410    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     411    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
     412    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
     413
     414    found  = .TRUE.
     415   
     416    SELECT CASE ( TRIM( variable ) )
     417!
     418!--    s grid
     419       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz'  )
     420
     421          grid_x = 'x'
     422          grid_y = 'y'
     423          grid_z = 'zu'   
     424!
     425!--    u grid
     426       CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
     427
     428          grid_x = 'xu'
     429          grid_y = 'y'
     430          grid_z = 'zu'
     431!
     432!--    v grid
     433       CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz'  )
     434
     435          grid_x = 'x'
     436          grid_y = 'yv'
     437          grid_z = 'zu'
     438!
     439!--    w grid
     440       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz'  )
     441
     442          grid_x = 'x'
     443          grid_y = 'y'
     444          grid_z = 'zw'
     445
     446       CASE DEFAULT
     447          found  = .FALSE.
     448          grid_x = 'none'
     449          grid_y = 'none'
     450          grid_z = 'none'
     451
     452    END SELECT
     453
     454
     455 END SUBROUTINE doq_define_netcdf_grid
     456 
     457!------------------------------------------------------------------------------!
     458!
     459! Description:
     460! ------------
     461!> Subroutine defining 2D output variables
     462!------------------------------------------------------------------------------!
     463 SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
     464                           mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
     465
     466
     467    IMPLICIT NONE
     468
     469    CHARACTER (LEN=*) ::  grid     !<
     470    CHARACTER (LEN=*) ::  mode     !<
     471    CHARACTER (LEN=*) ::  variable !<
     472
     473    INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
     474    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
     475    INTEGER(iwp) ::  i        !< grid index x-direction
     476    INTEGER(iwp) ::  j        !< grid index y-direction
     477    INTEGER(iwp) ::  k        !< grid index z-direction
     478    INTEGER(iwp) ::  nzb_do   !<
     479    INTEGER(iwp) ::  nzt_do   !<
     480
     481    LOGICAL ::  found             !< true if variable is in list
     482    LOGICAL ::  resorted          !< true if array is resorted
     483    LOGICAL ::  two_d             !< flag parameter that indicates 2D variables (horizontal cross sections)
     484
     485    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
     486   
     487    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
     488    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
     489   
     490    flag_nr  = 0
     491    found    = .TRUE.
     492    resorted = .FALSE.
     493    two_d    = .FALSE.
     494
     495    SELECT CASE ( TRIM( variable ) )
     496
     497       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
     498           IF ( av == 0 )  THEN
     499              to_be_resorted => ti
     500           ELSE
     501              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
     502                 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     503                 ti_av = REAL( fill_value, KIND = wp )
     504              ENDIF
     505              to_be_resorted => ti_av
     506           ENDIF
     507           flag_nr = 0
     508           
     509           IF ( mode == 'xy' )  grid = 'zu'
     510   
     511       CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
     512          IF ( av == 0 )  THEN
     513             to_be_resorted => uu
     514          ELSE
     515             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
     516                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     517                uu_av = REAL( fill_value, KIND = wp )
     518             ENDIF
     519             to_be_resorted => uu_av
     520          ENDIF
     521          flag_nr = 1
     522         
     523          IF ( mode == 'xy' )  grid = 'zu'
     524         
     525       CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
     526          IF ( av == 0 )  THEN
     527             to_be_resorted => vv
     528          ELSE
     529             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
     530                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     531                vv_av = REAL( fill_value, KIND = wp )
     532             ENDIF
     533             to_be_resorted => vv_av
     534          ENDIF
     535          flag_nr = 2
     536         
     537          IF ( mode == 'xy' )  grid = 'zu'
     538               
     539       CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
     540          IF ( av == 0 )  THEN
     541             to_be_resorted => ww
     542          ELSE
     543             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
     544                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     545                ww_av = REAL( fill_value, KIND = wp )
     546             ENDIF
     547             to_be_resorted => ww_av
     548          ENDIF
     549          flag_nr = 3
     550         
     551          IF ( mode == 'xy' )  grid = 'zw'
     552
     553       CASE DEFAULT
     554          found = .FALSE.
     555          grid  = 'none'
     556
     557    END SELECT
     558   
     559    IF ( found  .AND.  .NOT. resorted )  THEN     
     560       DO  i = nxl, nxr
     561          DO  j = nys, nyn
     562             DO  k = nzb_do, nzt_do
     563                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
     564                                         REAL( fill_value, KIND = wp ),        &
     565                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
     566             ENDDO
     567          ENDDO
     568       ENDDO
     569    ENDIF
     570 
     571 END SUBROUTINE doq_output_2d
     572 
     573 
     574!------------------------------------------------------------------------------!
     575!
     576! Description:
     577! ------------
     578!> Subroutine defining 3D output variables
     579!------------------------------------------------------------------------------!
     580 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
     581                           nzt_do )
     582 
     583    IMPLICIT NONE
     584
     585    CHARACTER (LEN=*) ::  variable !<
     586
     587    INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
     588    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
     589    INTEGER(iwp) ::  i        !< index variable along x-direction
     590    INTEGER(iwp) ::  j        !< index variable along y-direction
     591    INTEGER(iwp) ::  k        !< index variable along z-direction
     592    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
     593    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
     594
     595    LOGICAL ::  found             !< true if variable is in list
     596    LOGICAL ::  resorted          !< true if array is resorted
     597
     598    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
     599
     600    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !<
     601    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
     602
     603    flag_nr  = 0
     604    found    = .TRUE.
     605    resorted = .FALSE.
     606   
     607    SELECT CASE ( TRIM( variable ) )
     608
     609       CASE ( 'ti' )
     610          IF ( av == 0 )  THEN
     611             to_be_resorted => ti
     612          ELSE
     613             IF ( .NOT. ALLOCATED( ti_av ) ) THEN
     614                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     615                ti_av = REAL( fill_value, KIND = wp )
     616             ENDIF
     617             to_be_resorted => ti_av
     618          ENDIF
     619          flag_nr = 0
     620   
     621       CASE ( 'uu' )
     622          IF ( av == 0 )  THEN
     623             to_be_resorted => uu
     624          ELSE
     625             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
     626                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     627                uu_av = REAL( fill_value, KIND = wp )
     628             ENDIF
     629             to_be_resorted => uu_av
     630          ENDIF
     631          flag_nr = 1
     632             
     633       CASE ( 'vv' )
     634          IF ( av == 0 )  THEN
     635             to_be_resorted => vv
     636          ELSE
     637             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
     638                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     639                vv_av = REAL( fill_value, KIND = wp )
     640             ENDIF
     641             to_be_resorted => vv_av
     642          ENDIF
     643          flag_nr = 2
     644             
     645       CASE ( 'ww' )
     646          IF ( av == 0 )  THEN
     647             to_be_resorted => ww
     648          ELSE
     649             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
     650                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     651                ww_av = REAL( fill_value, KIND = wp )
     652             ENDIF
     653             to_be_resorted => ww_av
     654          ENDIF
     655          flag_nr = 3
     656
     657       CASE DEFAULT
     658          found = .FALSE.
     659
     660    END SELECT
     661   
     662    IF ( found  .AND.  .NOT. resorted )  THEN     
     663       DO  i = nxl, nxr
     664          DO  j = nys, nyn
     665             DO  k = nzb_do, nzt_do
     666                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
     667                                         REAL( fill_value, KIND = wp ),        &
     668                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
     669             ENDDO
     670          ENDDO
     671       ENDDO
     672    ENDIF
     673
     674 END SUBROUTINE doq_output_3d
     675 
     676! Description:
     677! ------------
     678!> Resorts the user-defined output quantity with indices (k,j,i) to a
     679!> temporary array with indices (i,j,k) for masked data output.
     680!------------------------------------------------------------------------------!
     681 SUBROUTINE doq_output_mask( av, variable, found, local_pf )
     682 
     683    USE control_parameters
     684       
     685    USE indices
     686   
     687    USE surface_mod,                                                           &
     688        ONLY:  get_topography_top_index_ji
     689 
     690    IMPLICIT NONE
     691
     692    CHARACTER (LEN=*) ::  variable  !<
     693    CHARACTER (LEN=5) ::  grid      !< flag to distinquish between staggered grids
     694
     695    INTEGER(iwp) ::  av           !< index indicating averaged or instantaneous output
     696    INTEGER(iwp) ::  flag_nr      !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
     697    INTEGER(iwp) ::  i            !< index variable along x-direction
     698    INTEGER(iwp) ::  j            !< index variable along y-direction
     699    INTEGER(iwp) ::  k            !< index variable along z-direction
     700    INTEGER(iwp) ::  topo_top_ind !< k index of highest horizontal surface
     701
     702    LOGICAL ::  found             !< true if variable is in list
     703    LOGICAL ::  resorted          !< true if array is resorted
     704
     705    REAL(wp),                                                                  &
     706       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
     707          local_pf   !<
     708    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
     709
     710    flag_nr  = 0
     711    found    = .TRUE.
     712    resorted = .FALSE.
     713    grid     = 's'
     714
     715    SELECT CASE ( TRIM( variable ) )
     716
     717       CASE ( 'ti' )
     718          IF ( av == 0 )  THEN
     719             to_be_resorted => ti
     720          ELSE
     721             to_be_resorted => ti_av
     722          ENDIF
     723          grid = 's'
     724          flag_nr = 0
     725   
     726       CASE ( 'uu' )
     727          IF ( av == 0 )  THEN
     728             to_be_resorted => uu
     729          ELSE
     730             to_be_resorted => uu_av
     731          ENDIF
     732          grid = 'u'
     733          flag_nr = 1
     734   
     735       CASE ( 'vv' )
     736          IF ( av == 0 )  THEN
     737             to_be_resorted => vv
     738          ELSE
     739             to_be_resorted => vv_av
     740          ENDIF
     741          grid = 'v'
     742          flag_nr = 2
     743   
     744       CASE ( 'ww' )
     745          IF ( av == 0 )  THEN
     746             to_be_resorted => ww
     747          ELSE
     748             to_be_resorted => ww_av
     749          ENDIF
     750          grid = 'w'
     751          flag_nr = 3
     752
     753
     754       CASE DEFAULT
     755          found = .FALSE.
     756
     757    END SELECT
     758   
     759    IF ( found  .AND.  .NOT. resorted )  THEN
     760       IF ( .NOT. mask_surface(mid) )  THEN
     761!
     762!--       Default masked output
     763          DO  i = 1, mask_size_l(mid,1)
     764             DO  j = 1, mask_size_l(mid,2)
     765                DO  k = 1, mask_size_l(mid,3)
     766                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k),            &
     767                                                     mask_j(mid,j),            &
     768                                                     mask_i(mid,i))
     769                ENDDO
     770             ENDDO
     771          ENDDO
     772
     773       ELSE
     774!
     775!--       Terrain-following masked output
     776          DO  i = 1, mask_size_l(mid,1)
     777             DO  j = 1, mask_size_l(mid,2)
     778!
     779!--             Get k index of highest horizontal surface
     780                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),     &
     781                                                            mask_i(mid,i),     &
     782                                                            grid )
     783!
     784!--             Save output array
     785                DO  k = 1, mask_size_l(mid,3)
     786                   local_pf(i,j,k) = to_be_resorted(                           &
     787                                          MIN( topo_top_ind+mask_k(mid,k),     &
     788                                               nzt+1 ),                        &
     789                                          mask_j(mid,j),                       &
     790                                          mask_i(mid,i)                     )
     791                ENDDO
     792             ENDDO
     793          ENDDO
     794
     795       ENDIF
     796    ENDIF
     797   
     798 END SUBROUTINE doq_output_mask
     799
     800!------------------------------------------------------------------------------!
     801! Description:
     802! ------------
     803!> Allocate required arrays
     804!------------------------------------------------------------------------------!
     805 SUBROUTINE doq_init
     806
     807    IMPLICIT NONE
     808   
     809    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
    158810!
    159811!-- Next line is to avoid compiler warnings about unused variables
     
    161813
    162814    initialized_diagnostic_output_quantities = .FALSE.
    163 
    164     IF ( .NOT. ALLOCATED( ti ) )  THEN
    165        ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
    166        ti = 0.0_wp
    167     ENDIF
     815   
     816    ivar = 1
     817   
     818    DO  WHILE ( ivar <= SIZE( do_all ) )
     819
     820       SELECT CASE ( TRIM( do_all(ivar) ) )
     821!
     822!--       Allocate array for 'turbulence intensity'
     823          CASE ( 'ti' )
     824             IF ( .NOT. ALLOCATED( ti ) )  THEN
     825                ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
     826                ti = 0.0_wp
     827             ENDIF
     828!
     829!--       Allocate array for uu
     830          CASE ( 'uu' )
     831             IF ( .NOT. ALLOCATED( uu ) )  THEN
     832                ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
     833                uu = 0.0_wp
     834             ENDIF
     835
     836!
     837!--       Allocate array for vv
     838          CASE ( 'vv' )
     839             IF ( .NOT. ALLOCATED( vv ) )  THEN
     840                ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
     841                vv = 0.0_wp
     842             ENDIF
     843
     844!
     845!--       Allocate array for ww
     846          CASE ( 'ww' )
     847             IF ( .NOT. ALLOCATED( ww ) )  THEN
     848                ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
     849                ww = 0.0_wp
     850             ENDIF
     851
     852       END SELECT
     853
     854       ivar = ivar + 1
     855    ENDDO
    168856
    169857    initialized_diagnostic_output_quantities = .TRUE.
    170 
    171  END SUBROUTINE diagnostic_output_quantities_init
     858   
     859 END SUBROUTINE doq_init
    172860
    173861
     
    175863! Description:
    176864! ------------
    177 !> ...
     865!> Calculate diagnostic quantities
    178866!--------------------------------------------------------------------------------------------------!
    179  SUBROUTINE diagnostic_output_quantities_calculate
    180 
    181 
    182     USE arrays_3d,                                                                                 &
    183         ONLY:  ddzu, u, v, w
    184 
    185     USE grid_variables,                                                                            &
    186         ONLY:  ddx, ddy
    187 
    188     USE indices,                                                                                   &
    189         ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
     867 SUBROUTINE doq_calculate
    190868
    191869    IMPLICIT NONE
    192870
    193871    INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
    194     INTEGER(iwp) ::  ivar_all   !< loop index over all 2d/3d/mask output quantities
     872    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
    195873
    196874
     
    198876!
    199877!-- Preparatory steps and initialization of output arrays
    200     IF ( .NOT.  prepared_diagnostic_output_quantities )     CALL diagnostic_output_quantities_prepare
    201     IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL diagnostic_output_quantities_init
    202 !
    203 !-- Save timestep number to check in time_integration if diagnostic_output_quantities_calculate
     878    IF ( .NOT.  prepared_diagnostic_output_quantities )     CALL doq_prepare
     879    IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL doq_init
     880!
     881!-- Save timestep number to check in time_integration if doq_calculate
    204882!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
    205883!-- done only once per timestep.
    206884    timestep_number_at_prev_calc = current_timestep_number
    207885
    208 
    209     ivar_all = 1
    210 
    211     DO  WHILE ( do_all(ivar_all)(1:1) /= ' ' )
    212 
    213        SELECT CASE ( TRIM( do_all(ivar_all) ) )
     886    ivar = 1
     887
     888    DO  WHILE ( ivar <= SIZE( do_all ) )
     889
     890       SELECT CASE ( TRIM( do_all(ivar) ) )
    214891!
    215892!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
     
    218895                DO  j = nys, nyn
    219896                   DO  k = nzb+1, nzt
    220 
    221                       ti(k,j,i) = 0.25_wp * SQRT(                                                  &
    222                         (   (   w(k,j+1,i) + w(k-1,j+1,i)                                          &
    223                               - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                                  &
    224                           - (   v(k+1,j,i) + v(k+1,j+1,i)                                          &
    225                               - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2                         &
    226                       + (   (   u(k+1,j,i) + u(k+1,j,i+1)                                          &
    227                               - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)                              &
    228                           - (   w(k,j,i+1) + w(k-1,j,i+1)                                          &
    229                               - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2                         &
    230                       + (   (   v(k,j,i+1) + v(k,j+1,i+1)                                          &
    231                               - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                                  &
    232                           - (   u(k,j+1,i) + u(k,j+1,i+1)                                          &
    233                               - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )                      &
    234                                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
     897                      ti(k,j,i) = 0.25_wp * SQRT(                              &
     898                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
     899                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
     900                          - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
     901                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
     902                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
     903                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
     904                          - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
     905                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
     906                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
     907                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
     908                          - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
     909                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
     910                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
    235911                   ENDDO
    236912                ENDDO
     913             ENDDO           
     914!
     915!--       uu
     916          CASE ( 'uu' )
     917             DO  i = nxl, nxr
     918                DO  j = nys, nyn
     919                   DO  k = nzb+1, nzt
     920                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
     921                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) )
     922                   ENDDO
     923                ENDDO
    237924             ENDDO
    238 
    239           END SELECT
    240 
    241           ivar_all = ivar_all + 1
     925!
     926!--       vv
     927          CASE ( 'vv' )
     928             DO  i = nxl, nxr
     929                DO  j = nys, nyn
     930                   DO  k = nzb+1, nzt
     931                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
     932                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) )
     933                   ENDDO
     934                ENDDO
     935             ENDDO
     936!
     937!--       ww
     938          CASE ( 'ww' )
     939             DO  i = nxl, nxr
     940                DO  j = nys, nyn
     941                   DO  k = nzb+1, nzt-1
     942                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
     943                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) )
     944                   ENDDO
     945                ENDDO
     946             ENDDO
     947
     948       END SELECT
     949
     950       ivar = ivar + 1
    242951    ENDDO
    243952
    244953!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
    245954
    246  END SUBROUTINE diagnostic_output_quantities_calculate
     955 END SUBROUTINE doq_calculate
    247956
    248957
     
    252961!> ...
    253962!------------------------------------------------------------------------------!
    254  SUBROUTINE diagnostic_output_quantities_prepare
    255 
    256 
    257     USE control_parameters,                                                                        &
     963 SUBROUTINE doq_prepare
     964
     965
     966    USE control_parameters,                                                    &
    258967        ONLY:  do2d, do3d, domask, masks, mid
    259968
     
    268977    INTEGER(iwp) ::  l          !< index for cutting string
    269978
    270 
    271979    prepared_diagnostic_output_quantities = .FALSE.
    272980
     
    280988       do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
    281989!
    282 !--    Gather 2d output quantity names, check for double occurrence of output quantity
     990!--    Gather 2d output quantity names.
     991!--    Check for double occurrence of output quantity, e.g. by _xy,
     992!--    _yz, _xz.
    283993       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
    284994          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
     
    2931003       ivar = 1
    2941004!
    295 !--    Gather 3d output quantity names, check for double occurrence of output quantity
     1005!--    Gather 3d output quantity names
    2961006       DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
    2971007          do_all(ivar_all) = do3d(av,ivar)
    298 
    2991008          ivar = ivar + 1
    3001009          ivar_all = ivar_all + 1
     
    3031012       ivar = 1
    3041013!
    305 !--    Gather masked output quantity names, check for double occurrence of output quantity
     1014!--    Gather masked output quantity names. Also check for double output
     1015!--    e.g. by different masks.
    3061016       DO  mid = 1, masks
    3071017          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
    308              do_all(ivar_all) = domask(mid,av,ivar)
     1018             IF ( .NOT.  ANY( do_all == domask(mid,av,ivar) ) )  THEN
     1019                do_all(ivar_all) = do2d_var(av,ivar)
     1020             ENDIF
    3091021
    3101022             ivar = ivar + 1
     
    3181030    prepared_diagnostic_output_quantities = .TRUE.
    3191031
    320  END SUBROUTINE diagnostic_output_quantities_prepare
     1032 END SUBROUTINE doq_prepare
     1033 
     1034!------------------------------------------------------------------------------!
     1035! Description:
     1036! ------------
     1037!> Subroutine reads local (subdomain) restart data
     1038!> Note: With the current structure reading of non-standard array is not
     1039!> possible
     1040!------------------------------------------------------------------------------!
     1041!  SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
     1042!                            nxr_on_file, nynf, nync, nyn_on_file, nysf,         &
     1043!                            nysc, nys_on_file, tmp_3d_non_standard, found )
     1044
     1045!
     1046!     USE control_parameters
     1047!         
     1048!     USE indices
     1049!     
     1050!     USE kinds
     1051!     
     1052!     USE pegrid
     1053!
     1054!
     1055!     IMPLICIT NONE
     1056!
     1057!     INTEGER(iwp) ::  k               !<
     1058!     INTEGER(iwp) ::  nxlc            !<
     1059!     INTEGER(iwp) ::  nxlf            !<
     1060!     INTEGER(iwp) ::  nxl_on_file     !<
     1061!     INTEGER(iwp) ::  nxrc            !<
     1062!     INTEGER(iwp) ::  nxrf            !<
     1063!     INTEGER(iwp) ::  nxr_on_file     !<
     1064!     INTEGER(iwp) ::  nync            !<
     1065!     INTEGER(iwp) ::  nynf            !<
     1066!     INTEGER(iwp) ::  nyn_on_file     !<
     1067!     INTEGER(iwp) ::  nysc            !<
     1068!     INTEGER(iwp) ::  nysf            !<
     1069!     INTEGER(iwp) ::  nys_on_file     !<
     1070!
     1071!     LOGICAL, INTENT(OUT)  :: found
     1072!     
     1073!     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
     1074! !
     1075! !-- If temporary non-standard array for reading is already allocated,
     1076! !-- deallocate it.
     1077!     IF ( ALLOCATED( tmp_3d_non_standard ) )  DEALLOCATE( tmp_3d_non_standard )
     1078!     
     1079!     found = .TRUE.
     1080!
     1081!     SELECT CASE ( restart_string(1:length) )
     1082!
     1083!        CASE ( 'ti_av' )
     1084!           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
     1085!              ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     1086!           ENDIF
     1087!           IF ( k == 1 )  THEN
     1088!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
     1089!                                            nxl_on_file:nxr_on_file) )
     1090!              READ ( 13 )  tmp_3d_non_standard
     1091!           ENDIF
     1092!           ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
     1093!     
     1094!        CASE ( 'uu_av' )
     1095!           IF ( .NOT. ALLOCATED( uu_av ) )  THEN
     1096!              ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     1097!           ENDIF
     1098!           IF ( k == 1 )  THEN
     1099!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
     1100!                                            nxl_on_file:nxr_on_file) )
     1101!              READ ( 13 )  tmp_3d_non_standard
     1102!           ENDIF
     1103!           uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
     1104!                   
     1105!        CASE ( 'vv_av' )
     1106!           IF ( .NOT. ALLOCATED( vv_av ) )  THEN
     1107!              ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     1108!           ENDIF
     1109!           IF ( k == 1 )  THEN
     1110!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
     1111!                                            nxl_on_file:nxr_on_file) )
     1112!              READ ( 13 )  tmp_3d_non_standard
     1113!           ENDIF
     1114!           vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
     1115!                   
     1116!        CASE ( 'ww_av' )
     1117!           IF ( .NOT. ALLOCATED( ww_av ) )  THEN
     1118!              ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     1119!           ENDIF
     1120!           IF ( k == 1 )  THEN
     1121!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
     1122!                                            nxl_on_file:nxr_on_file) )
     1123!              READ ( 13 )  tmp_3d_non_standard
     1124!           ENDIF
     1125!           ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
     1126!                         
     1127!
     1128!        CASE DEFAULT
     1129!
     1130!           found = .FALSE.
     1131!
     1132!     END SELECT
     1133!
     1134!  END SUBROUTINE doq_rrd_local
     1135 
     1136!------------------------------------------------------------------------------!
     1137! Description:
     1138! ------------
     1139!> Subroutine writes local (subdomain) restart data
     1140!------------------------------------------------------------------------------!
     1141 SUBROUTINE doq_wrd_local
     1142
     1143
     1144    IMPLICIT NONE
     1145
     1146    IF ( ALLOCATED( ti_av ) )  THEN
     1147       CALL wrd_write_string( 'ti_av' )
     1148       WRITE ( 14 )  ti_av
     1149    ENDIF
     1150   
     1151    IF ( ALLOCATED( uu_av ) )  THEN
     1152       CALL wrd_write_string( 'uu_av' )
     1153       WRITE ( 14 )  uu_av
     1154    ENDIF
     1155       
     1156    IF ( ALLOCATED( vv_av ) )  THEN
     1157       CALL wrd_write_string( 'vv_av' )
     1158       WRITE ( 14 )  vv_av
     1159    ENDIF
     1160       
     1161    IF ( ALLOCATED( ww_av ) )  THEN
     1162       CALL wrd_write_string( 'ww_av' )
     1163       WRITE ( 14 )  ww_av
     1164    ENDIF
     1165
     1166
     1167 END SUBROUTINE doq_wrd_local
     1168 
     1169 
    3211170
    3221171 END MODULE diagnostic_output_quantities_mod
Note: See TracChangeset for help on using the changeset viewer.