Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 4583 for palm


Ignore:
Timestamp:
Jun 29, 2020 12:36:47 PM (4 years ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/data_output_tseries.f90

    r4360 r4583  
    11!> @file data_output_tseries.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    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.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    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/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
    2729! Corrected "Former revisions" section
    28 ! 
     30!
    2931! 3655 2019-01-07 16:51:22Z knoop
    3032! unused format removed
    … …  
    3638! Description:
    3739! ------------
    38 !> Time series output for PROFIL. Always all time series are stored. A selection
    39 !> can be applied via the PROFIL-parameters in close_file.
    40 !------------------------------------------------------------------------------!
     40!> Time series output for PROFIL. Always all time series are stored. A selection can be applied via
     41!> the PROFIL-parameters in close_file.
     42!--------------------------------------------------------------------------------------------------!
    4143 SUBROUTINE data_output_tseries
    42  
    4344
    44     USE control_parameters,                                                    &
     45
     46    USE control_parameters,                                                                        &
    4547        ONLY:  dots_time_count, time_since_reference_point
    4648
    47     USE cpulog,                                                                &
    48         ONLY:  cpu_log, log_point 
     49    USE cpulog,                                                                                    &
     50        ONLY:  cpu_log, log_point
    4951
    5052    USE kinds
    … …  
    5355    USE NETCDF
    5456#endif
    55     USE netcdf_interface,                                                      &
    56         ONLY:  dots_num, id_set_ts, id_var_dots, id_var_time_ts, nc_stat,      &
    57                netcdf_handle_error
     57    USE netcdf_interface,                                                                          &
     58        ONLY:  dots_num, id_set_ts, id_var_dots, id_var_time_ts, nc_stat, netcdf_handle_error
    5859
    5960    USE pegrid
    6061
    6162    USE profil_parameter
    62    
    63     USE statistics,                                                            &
     63
     64    USE statistics,                                                                                &
    6465        ONLY:  flow_statistics_called, statistic_regions, ts_value
    6566
    … …  
    8485!--    Open file for time series output in NetCDF format
    8586       CALL check_open( 105 )
    86        
     87
    8788!--    Increment the counter for number of output times
    88 !      CAUTION: The following line has to be after the call of the subroutine
    89 !               check_open, since check_open resets the counter dots_time_count
    90 !               to 0, if a new file is opened
     89!--    CAUTION: The following line has to be after the call of the subroutine check_open, since
     90!--             check_open resets the counter dots_time_count to 0, if a new file is opened
    9191       dots_time_count = dots_time_count + 1
    92        
     92
    9393#if defined( __netcdf )
    9494!
    9595!--    Update the time series time axis
    96        nc_stat = NF90_PUT_VAR( id_set_ts, id_var_time_ts,        &
    97                                (/ time_since_reference_point /), &
    98                                start = (/ dots_time_count /),    &
     96       nc_stat = NF90_PUT_VAR( id_set_ts, id_var_time_ts,                                          &
     97                               (/ time_since_reference_point /),                                   &
     98                               start = (/ dots_time_count /),                                      &
    9999                               count = (/ 1 /) )
    100100       CALL netcdf_handle_error( 'data_output_tseries', 350 )
    … …  
    108108#if defined( __netcdf )
    109109          DO  i = 1, dots_num
    110              nc_stat = NF90_PUT_VAR( id_set_ts, id_var_dots(i,sr),  &
    111                                      (/ ts_value(i,sr) /),          &
    112                                      start = (/ dots_time_count /), &
     110             nc_stat = NF90_PUT_VAR( id_set_ts, id_var_dots(i,sr),                                 &
     111                                     (/ ts_value(i,sr) /),                                         &
     112                                     start = (/ dots_time_count /),                                &
    113113                                     count = (/ 1 /) )
    114114             CALL netcdf_handle_error( 'data_output_tseries', 351 )

  • palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90

    r4560 r4583  
    11!> @file diagnostic_output_quantities_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4560 2020-06-11 12:19:47Z suehring
    2729! - Bugfix in calculation of vertical momentum and scalar fluxes
    2830! - remove averaged output variables from PUBLIC list
    29 ! 
     31!
    3032! 4535 2020-05-15 12:07:23Z raasch
    3133! bugfix for restart data format query
    32 ! 
     34!
    3335! 4518 2020-05-04 15:44:28Z suehring
    34 ! * Define arrays over ghost points in order to allow for standard mpi-io
    35 !   treatment. By this modularization of restart-data input is possible with
    36 !   the module interface.
     36! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this
     37!   modularization of restart-data input is possible with the module interface.
    3738! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.
    38 ! 
     39!
    3940! 4517 2020-05-03 14:29:30Z raasch
    4041! use statement for exchange horiz added,
    4142! bugfix for call of exchange horiz 2d
    42 ! 
     43!
    4344! 4431 2020-02-27 23:23:01Z gronemeier
    4445! added wspeed and wdir output; bugfix: set fill_value in case of masked output
    4546!
    4647! 4360 2020-01-07 11:25:50Z suehring
    47 ! added output of wu, wv, wtheta and wq to enable covariance calculation
    48 ! according to temporal EC method
     48! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC
     49! method
    4950!
    5051! 4346 2019-12-18 11:55:56Z motisi
    51 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    52 ! topography information used in wall_flags_static_0
     52! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     53! information used in wall_flags_static_0
    5354!
    5455! 4331 2019-12-10 18:25:02Z suehring
    … …  
    6364!
    6465! 4167 2019-08-16 11:01:48Z suehring
    65 ! Changed behaviour of masked output over surface to follow terrain and ignore
    66 ! buildings (J.Resler, T.Gronemeier)
     66! Changed behaviour of masked output over surface to follow terrain and ignore buildings
     67! (J.Resler, T.Gronemeier)
    6768!
    6869! 4157 2019-08-14 09:19:12Z suehring
    69 ! Initialization restructured, in order to work also when data output during
    70 ! spin-up is enabled.
     70! Initialization restructured, in order to work also when data output during spin-up is enabled.
    7171!
    7272! 4132 2019-08-02 12:34:17Z suehring
    … …  
    7474!
    7575! 4069 2019-07-01 14:05:51Z Giersch
    76 ! Masked output running index mid has been introduced as a local variable to
    77 ! avoid runtime error (Loop variable has been modified) in time_integration
     76! Masked output running index mid has been introduced as a local variable to avoid runtime error
     77! (Loop variable has been modified) in time_integration
    7878!
    7979! 4039 2019-06-18 10:32:41Z suehring
    80 ! - Add output of uu, vv, ww to enable variance calculation according temporal
    81 !   EC method
     80! - Add output of uu, vv, ww to enable variance calculation according temporal EC method
    8281! - Allocate arrays only when they are required
    8382! - Formatting adjustment
    … …  
    8988!
    9089! 3995 2019-05-22 18:59:54Z suehring
    91 ! Avoid compiler warnings about unused variable and fix string operation which
    92 ! is not allowed with PGI compiler
     90! Avoid compiler warnings about unused variable and fix string operation which is not allowed with
     91! PGI compiler
    9392!
    9493! 3994 2019-05-22 18:08:09Z suehring
    … …  
    103102! ------------
    104103!> ...
    105 !------------------------------------------------------------------------------!
     104!--------------------------------------------------------------------------------------------------!
    106105 MODULE diagnostic_output_quantities_mod
    107106
    108     USE arrays_3d,                                                             &
    109         ONLY:  ddzu,                                                           &
    110                pt,                                                             &
    111                q,                                                              &
    112                u,                                                              &
    113                v,                                                              &
    114                w,                                                              &
    115                zu,                                                             &
     107    USE arrays_3d,                                                                                 &
     108        ONLY:  ddzu,                                                                               &
     109               pt,                                                                                 &
     110               q,                                                                                  &
     111               u,                                                                                  &
     112               v,                                                                                  &
     113               w,                                                                                  &
     114               zu,                                                                                 &
    116115               zw
    117116
    118     USE basic_constants_and_equations_mod,                                     &
     117    USE basic_constants_and_equations_mod,                                                         &
    119118        ONLY:  kappa, pi
    120119
    121     USE control_parameters,                                                    &
    122         ONLY:  current_timestep_number,                                        &
    123                data_output,                                                    &
    124                message_string,                                                 &
    125                restart_data_format_output,                                     &
     120    USE control_parameters,                                                                        &
     121        ONLY:  current_timestep_number,                                                            &
     122               data_output,                                                                        &
     123               message_string,                                                                     &
     124               restart_data_format_output,                                                         &
    126125               varnamelength
    127126!
    128 !     USE cpulog,                                                                &
     127!     USE cpulog,                                                                                   &
    129128!         ONLY:  cpu_log, log_point
    130129
    131     USE exchange_horiz_mod,                                                    &
     130    USE exchange_horiz_mod,                                                                        &
    132131        ONLY:  exchange_horiz_2d
    133132
    134    USE grid_variables,                                                         &
     133   USE grid_variables,                                                                             &
    135134        ONLY:  ddx, ddy
    136135
    137     USE indices,                                                               &
    138         ONLY:  nbgp,                                                           &
    139                nxl,                                                            &
    140                nxlg,                                                           &
    141                nxr,                                                            &
    142                nxrg,                                                           &
    143                nyn,                                                            &
    144                nyng,                                                           &
    145                nys,                                                            &
    146                nysg,                                                           &
    147                nzb,                                                            &
    148                nzt,                                                            &
     136    USE indices,                                                                                   &
     137        ONLY:  nbgp,                                                                               &
     138               nxl,                                                                                &
     139               nxlg,                                                                               &
     140               nxr,                                                                                &
     141               nxrg,                                                                               &
     142               nyn,                                                                                &
     143               nyng,                                                                               &
     144               nys,                                                                                &
     145               nysg,                                                                               &
     146               nzb,                                                                                &
     147               nzt,                                                                                &
    149148               wall_flags_total_0
    150149
    … …  
    156155               wrd_mpi_io
    157156
    158     USE surface_mod,                                                           &
    159         ONLY:  surf_def_h,                                                     &
    160                surf_lsm_h,                                                     &
    161                surf_type,                                                      &
     157    USE surface_mod,                                                                               &
     158        ONLY:  surf_def_h,                                                                         &
     159               surf_lsm_h,                                                                         &
     160               surf_type,                                                                          &
    162161               surf_usm_h
    163162
    … …  
    209208!
    210209!-- Public variables
    211     PUBLIC do_all,                                                             &
    212            initialized_diagnostic_output_quantities,                           &
    213            prepared_diagnostic_output_quantities,                              &
     210    PUBLIC do_all,                                                                                 &
     211           initialized_diagnostic_output_quantities,                                               &
     212           prepared_diagnostic_output_quantities,                                                  &
    214213           timestep_number_at_prev_calc
    215214!
    216215!-- Public routines
    217     PUBLIC doq_3d_data_averaging,                                              &
    218            doq_calculate,                                                      &
    219            doq_check_data_output,                                              &
    220            doq_define_netcdf_grid,                                             &
    221            doq_init,                                                           &
    222            doq_output_2d,                                                      &
    223            doq_output_3d,                                                      &
    224            doq_output_mask,                                                    &
    225            doq_rrd_local,                                                      &
     216    PUBLIC doq_3d_data_averaging,                                                                  &
     217           doq_calculate,                                                                          &
     218           doq_check_data_output,                                                                  &
     219           doq_define_netcdf_grid,                                                                 &
     220           doq_init,                                                                               &
     221           doq_output_2d,                                                                          &
     222           doq_output_3d,                                                                          &
     223           doq_output_mask,                                                                        &
     224           doq_rrd_local,                                                                          &
    226225           doq_wrd_local
    227226
    … …  
    275274 CONTAINS
    276275
    277 !------------------------------------------------------------------------------!
     276!--------------------------------------------------------------------------------------------------!
    278277! Description:
    279278! ------------
    280 !> Sum up and time-average diagnostic output quantities as well as allocate
    281 !> the array necessary for storing the average.
    282 !------------------------------------------------------------------------------!
     279!> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for
     280!> storing the average.
     281!--------------------------------------------------------------------------------------------------!
    283282 SUBROUTINE doq_3d_data_averaging( mode, variable )
    284283
    285     USE control_parameters,                                                    &
     284    USE control_parameters,                                                                        &
    286285        ONLY:  average_count_3d
    287286
    … …  
    293292    INTEGER(iwp) ::  k !<
    294293
     294
    295295    IF ( mode == 'allocate' )  THEN
    296296
    … …  
    383383
    384384          CASE ( 'ti' )
    385              IF ( ALLOCATED( ti_av ) ) THEN
     385             IF ( ALLOCATED( ti_av ) )  THEN
    386386                DO  i = nxl, nxr
    387387                   DO  j = nys, nyn
    … …  
    394394
    395395          CASE ( 'uu' )
    396              IF ( ALLOCATED( uu_av ) ) THEN
     396             IF ( ALLOCATED( uu_av ) )  THEN
    397397                DO  i = nxl, nxr
    398398                   DO  j = nys, nyn
    … …  
    405405
    406406          CASE ( 'vv' )
    407              IF ( ALLOCATED( vv_av ) ) THEN
     407             IF ( ALLOCATED( vv_av ) )  THEN
    408408                DO  i = nxl, nxr
    409409                   DO  j = nys, nyn
    … …  
    416416
    417417          CASE ( 'ww' )
    418              IF ( ALLOCATED( ww_av ) ) THEN
     418             IF ( ALLOCATED( ww_av ) )  THEN
    419419                DO  i = nxl, nxr
    420420                   DO  j = nys, nyn
    … …  
    427427
    428428          CASE ( 'wu' )
    429              IF ( ALLOCATED( wu_av ) ) THEN
     429             IF ( ALLOCATED( wu_av ) )  THEN
    430430                DO  i = nxl, nxr
    431431                   DO  j = nys, nyn
    … …  
    438438
    439439          CASE ( 'wv' )
    440              IF ( ALLOCATED( wv_av ) ) THEN
     440             IF ( ALLOCATED( wv_av ) )  THEN
    441441                DO  i = nxl, nxr
    442442                   DO  j = nys, nyn
    … …  
    449449
    450450          CASE ( 'wtheta' )
    451              IF ( ALLOCATED( wtheta_av ) ) THEN
     451             IF ( ALLOCATED( wtheta_av ) )  THEN
    452452                DO  i = nxl, nxr
    453453                   DO  j = nys, nyn
    … …  
    460460
    461461          CASE ( 'wq' )
    462              IF ( ALLOCATED( wq_av ) ) THEN
     462             IF ( ALLOCATED( wq_av ) )  THEN
    463463                DO  i = nxl, nxr
    464464                   DO  j = nys, nyn
    … …  
    471471
    472472          CASE ( 'theta_2m*' )
    473              IF ( ALLOCATED( pt_2m_av ) ) THEN
     473             IF ( ALLOCATED( pt_2m_av ) )  THEN
    474474                DO  i = nxl, nxr
    475475                   DO  j = nys, nyn
    … …  
    480480
    481481          CASE ( 'wspeed_10m*' )
    482              IF ( ALLOCATED( uv_10m_av ) ) THEN
     482             IF ( ALLOCATED( uv_10m_av ) )  THEN
    483483                DO  i = nxl, nxr
    484484                   DO  j = nys, nyn
    … …  
    489489
    490490          CASE ( 'wspeed' )
    491             IF ( ALLOCATED( wspeed_av ) ) THEN
     491            IF ( ALLOCATED( wspeed_av ) )  THEN
    492492               DO  i = nxl, nxr
    493493                  DO  j = nys, nyn
    … …  
    500500
    501501          CASE ( 'wdir' )
    502              IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
     502             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
    503503                DO  i = nxl, nxr
    504504                   DO  j = nys, nyn
    … …  
    521521
    522522          CASE ( 'ti' )
    523              IF ( ALLOCATED( ti_av ) ) THEN
     523             IF ( ALLOCATED( ti_av ) )  THEN
    524524                DO  i = nxl, nxr
    525525                   DO  j = nys, nyn
    … …  
    532532
    533533          CASE ( 'uu' )
    534              IF ( ALLOCATED( uu_av ) ) THEN
     534             IF ( ALLOCATED( uu_av ) )  THEN
    535535                DO  i = nxl, nxr
    536536                   DO  j = nys, nyn
    … …  
    543543
    544544          CASE ( 'vv' )
    545              IF ( ALLOCATED( vv_av ) ) THEN
     545             IF ( ALLOCATED( vv_av ) )  THEN
    546546                DO  i = nxl, nxr
    547547                   DO  j = nys, nyn
    … …  
    554554
    555555          CASE ( 'ww' )
    556              IF ( ALLOCATED( ww_av ) ) THEN
     556             IF ( ALLOCATED( ww_av ) )  THEN
    557557                DO  i = nxl, nxr
    558558                   DO  j = nys, nyn
    … …  
    565565
    566566          CASE ( 'wu' )
    567              IF ( ALLOCATED( wu_av ) ) THEN
     567             IF ( ALLOCATED( wu_av ) )  THEN
    568568                DO  i = nxl, nxr
    569569                   DO  j = nys, nyn
    … …  
    576576
    577577          CASE ( 'wv' )
    578              IF ( ALLOCATED( wv_av ) ) THEN
     578             IF ( ALLOCATED( wv_av ) )  THEN
    579579                DO  i = nxl, nxr
    580580                   DO  j = nys, nyn
    … …  
    587587
    588588          CASE ( 'wtheta' )
    589              IF ( ALLOCATED( wtheta_av ) ) THEN
     589             IF ( ALLOCATED( wtheta_av ) )  THEN
    590590                DO  i = nxl, nxr
    591591                   DO  j = nys, nyn
    … …  
    598598
    599599          CASE ( 'wq' )
    600              IF ( ALLOCATED( wq_av ) ) THEN
     600             IF ( ALLOCATED( wq_av ) )  THEN
    601601                DO  i = nxl, nxr
    602602                   DO  j = nys, nyn
    … …  
    609609
    610610         CASE ( 'theta_2m*' )
    611             IF ( ALLOCATED( pt_2m_av ) ) THEN
     611            IF ( ALLOCATED( pt_2m_av ) )  THEN
    612612               DO  i = nxlg, nxrg
    613613                  DO  j = nysg, nyng
    … …  
    619619
    620620         CASE ( 'wspeed_10m*' )
    621             IF ( ALLOCATED( uv_10m_av ) ) THEN
     621            IF ( ALLOCATED( uv_10m_av ) )  THEN
    622622               DO  i = nxlg, nxrg
    623623                  DO  j = nysg, nyng
    … …  
    629629
    630630         CASE ( 'wspeed' )
    631              IF ( ALLOCATED( wspeed_av ) ) THEN
     631             IF ( ALLOCATED( wspeed_av ) )  THEN
    632632                DO  i = nxl, nxr
    633633                   DO  j = nys, nyn
    … …  
    640640
    641641          CASE ( 'wdir' )
    642              IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) ) THEN
     642             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
    643643
    644644                IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
    … …  
    652652                         u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    653653                         v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    654                          wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &
    655                                         / pi * 180.0_wp + 180.0_wp
     654                         wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) )          &
     655                                          / pi * 180.0_wp + 180.0_wp
    656656                      ENDDO
    657657                   ENDDO
    … …  
    666666 END SUBROUTINE doq_3d_data_averaging
    667667
    668 !------------------------------------------------------------------------------!
     668!--------------------------------------------------------------------------------------------------!
    669669! Description:
    670670! ------------
    671671!> Check data output for diagnostic output
    672 !------------------------------------------------------------------------------!
     672!--------------------------------------------------------------------------------------------------!
    673673 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
    674674
    … …  
    682682    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
    683683
     684
    684685    SELECT CASE ( TRIM( var ) )
    685686
    … …  
    719720!--       Check if output quantity is _xy only.
    720721          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    721              message_string = 'illegal value for data_output: "' //            &
    722                               TRIM( var ) // '" & only 2d-horizontal ' //      &
     722             message_string = 'illegal value for data_output: "' //                                &
     723                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
    723724                              'cross sections are allowed for this value'
    724725             CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
    … …  
    736737 END SUBROUTINE doq_check_data_output
    737738
    738 !------------------------------------------------------------------------------!
     739!--------------------------------------------------------------------------------------------------!
    739740!
    740741! Description:
    741742! ------------
    742743!> Subroutine defining appropriate grid for netcdf variables.
    743 !------------------------------------------------------------------------------!
     744!--------------------------------------------------------------------------------------------------!
    744745 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
    745746
    … …  
    747748
    748749    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
    749     LOGICAL, INTENT(OUT)           ::  found       !<
    750750    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
    751751    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
    752752    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
    753753
     754    LOGICAL, INTENT(OUT)           ::  found       !<
     755
     756
    754757    found  = .TRUE.
    755758
    … …  
    757760!
    758761!--    s grid
    759        CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                 &
    760               'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                 &
     762       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
     763              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
    761764              'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
    762765
    … …  
    787790!
    788791!--    w grid
    789        CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz',                                 &
    790               'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                 &
    791               'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                 &
    792               'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                 &
    793               'wq', 'wq_xy', 'wq_xz', 'wq_yz'  )
     792       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz',                                                     &
     793              'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                                     &
     794              'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                                     &
     795              'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                                     &
     796              'wq', 'wq_xy', 'wq_xz', 'wq_yz' )
    794797
    795798          grid_x = 'x'
    … …  
    808811 END SUBROUTINE doq_define_netcdf_grid
    809812
    810 !------------------------------------------------------------------------------!
     813!--------------------------------------------------------------------------------------------------!
    811814!
    812815! Description:
    813816! ------------
    814817!> Subroutine defining 2D output variables
    815 !------------------------------------------------------------------------------!
    816  SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
    817                            mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
     818!--------------------------------------------------------------------------------------------------!
     819 SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do,       &
     820                           fill_value )
    818821
    819822
    … …  
    841844    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
    842845
     846
    843847    flag_nr  = 0
    844848    found    = .TRUE.
    … …  
    852856              to_be_resorted => ti
    853857           ELSE
    854               IF ( .NOT. ALLOCATED( ti_av ) ) THEN
     858              IF ( .NOT. ALLOCATED( ti_av ) )  THEN
    855859                 ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    856860                 ti_av = REAL( fill_value, KIND = wp )
    … …  
    866870             to_be_resorted => uu
    867871          ELSE
    868              IF ( .NOT. ALLOCATED( uu_av ) ) THEN
     872             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
    869873                ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    870874                uu_av = REAL( fill_value, KIND = wp )
    … …  
    880884             to_be_resorted => vv
    881885          ELSE
    882              IF ( .NOT. ALLOCATED( vv_av ) ) THEN
     886             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
    883887                ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    884888                vv_av = REAL( fill_value, KIND = wp )
    … …  
    894898             to_be_resorted => ww
    895899          ELSE
    896              IF ( .NOT. ALLOCATED( ww_av ) ) THEN
     900             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
    897901                ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    898902                ww_av = REAL( fill_value, KIND = wp )
    … …  
    908912             to_be_resorted => wu
    909913          ELSE
    910              IF ( .NOT. ALLOCATED( wu_av ) ) THEN
     914             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
    911915                ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    912916                wu_av = REAL( fill_value, KIND = wp )
    … …  
    922926             to_be_resorted => wv
    923927          ELSE
    924              IF ( .NOT. ALLOCATED( wv_av ) ) THEN
     928             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
    925929                ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    926930                wv_av = REAL( fill_value, KIND = wp )
    … …  
    936940             to_be_resorted => wtheta
    937941          ELSE
    938              IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
     942             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
    939943                ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    940944                wtheta_av = REAL( fill_value, KIND = wp )
    … …  
    950954             to_be_resorted => wq
    951955          ELSE
    952              IF ( .NOT. ALLOCATED( wq_av ) ) THEN
     956             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
    953957                ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    954958                wq_av = REAL( fill_value, KIND = wp )
    … …  
    968972             ENDDO
    969973          ELSE
    970              IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
     974             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
    971975                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
    972976                pt_2m_av = REAL( fill_value, KIND = wp )
    … …  
    990994             ENDDO
    991995          ELSE
    992              IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
     996             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
    993997                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
    994998                uv_10m_av = REAL( fill_value, KIND = wp )
    … …  
    10081012             to_be_resorted => wspeed
    10091013          ELSE
    1010              IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
     1014             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
    10111015                ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10121016                wspeed_av = REAL( fill_value, KIND = wp )
    … …  
    10221026             to_be_resorted => wdir
    10231027          ELSE
    1024              IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
     1028             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
    10251029                ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10261030                wdir_av = REAL( fill_value, KIND = wp )
    … …  
    10421046          DO  j = nys, nyn
    10431047             DO  k = nzb_do, nzt_do
    1044                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
    1045                                      REAL( fill_value, KIND = wp ),            &
    1046                                      BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
     1048                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
     1049                                         REAL( fill_value, KIND = wp ),                            &
     1050                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
    10471051             ENDDO
    10481052          ENDDO
    … …  
    10531057
    10541058
    1055 !------------------------------------------------------------------------------!
     1059!--------------------------------------------------------------------------------------------------!
    10561060!
    10571061! Description:
    10581062! ------------
    10591063!> Subroutine defining 3D output variables
    1060 !------------------------------------------------------------------------------!
    1061  SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
    1062                            nzt_do )
     1064!--------------------------------------------------------------------------------------------------!
     1065 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    10631066
    10641067    IMPLICIT NONE
    … …  
    10821085    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
    10831086
     1087
    10841088    flag_nr  = 0
    10851089    found    = .TRUE.
    … …  
    10921096             to_be_resorted => ti
    10931097          ELSE
    1094              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
     1098             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
    10951099                ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10961100                ti_av = REAL( fill_value, KIND = wp )
    … …  
    11041108             to_be_resorted => uu
    11051109          ELSE
    1106              IF ( .NOT. ALLOCATED( uu_av ) ) THEN
     1110             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
    11071111                ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11081112                uu_av = REAL( fill_value, KIND = wp )
    … …  
    11161120             to_be_resorted => vv
    11171121          ELSE
    1118              IF ( .NOT. ALLOCATED( vv_av ) ) THEN
     1122             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
    11191123                ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11201124                vv_av = REAL( fill_value, KIND = wp )
    … …  
    11281132             to_be_resorted => ww
    11291133          ELSE
    1130              IF ( .NOT. ALLOCATED( ww_av ) ) THEN
     1134             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
    11311135                ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11321136                ww_av = REAL( fill_value, KIND = wp )
    … …  
    11401144             to_be_resorted => wu
    11411145          ELSE
    1142              IF ( .NOT. ALLOCATED( wu_av ) ) THEN
     1146             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
    11431147                ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11441148                wu_av = REAL( fill_value, KIND = wp )
    … …  
    11521156             to_be_resorted => wv
    11531157          ELSE
    1154              IF ( .NOT. ALLOCATED( wv_av ) ) THEN
     1158             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
    11551159                ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11561160                wv_av = REAL( fill_value, KIND = wp )
    … …  
    11641168             to_be_resorted => wtheta
    11651169          ELSE
    1166              IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
     1170             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
    11671171                ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11681172                wtheta_av = REAL( fill_value, KIND = wp )
    … …  
    11761180             to_be_resorted => wq
    11771181          ELSE
    1178              IF ( .NOT. ALLOCATED( wq_av ) ) THEN
     1182             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
    11791183                ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11801184                wq_av = REAL( fill_value, KIND = wp )
    … …  
    11881192             to_be_resorted => wspeed
    11891193          ELSE
    1190              IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
     1194             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
    11911195                ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11921196                wspeed_av = REAL( fill_value, KIND = wp )
    … …  
    12001204             to_be_resorted => wdir
    12011205          ELSE
    1202              IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
     1206             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
    12031207                ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    12041208                wdir_av = REAL( fill_value, KIND = wp )
    … …  
    12171221          DO  j = nys, nyn
    12181222             DO  k = nzb_do, nzt_do
    1219                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
    1220                                      REAL( fill_value, KIND = wp ),            &
    1221                                      BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
     1223                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
     1224                                         REAL( fill_value, KIND = wp ),                            &
     1225                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
    12221226             ENDDO
    12231227          ENDDO
    … …  
    12271231 END SUBROUTINE doq_output_3d
    12281232
     1233
     1234!--------------------------------------------------------------------------------------------------!
     1235!
    12291236! Description:
    12301237! ------------
    1231 !> Resorts the user-defined output quantity with indices (k,j,i) to a
    1232 !> temporary array with indices (i,j,k) for masked data output.
    1233 !------------------------------------------------------------------------------!
     1238!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
     1239!> (i,j,k) for masked data output.
     1240!--------------------------------------------------------------------------------------------------!
    12341241 SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
    12351242
    … …  
    12391246
    12401247    IMPLICIT NONE
     1248
     1249    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
    12411250
    12421251    CHARACTER (LEN=*) ::  variable   !<
    … …  
    12571266    LOGICAL      ::  resorted        !< true if array is resorted
    12581267
    1259     REAL(wp),                                                                  &
    1260        DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
    1261           local_pf   !<
    1262     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
    1263 
    1264     REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
     1268    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
     1269
     1270    REAL(wp), DIMENSION(:,:,:), POINTER  ::  to_be_resorted  !< points to array which needs to be resorted for output
     1271
    12651272
    12661273    flag_nr  = 0
    … …  
    13731380             DO  j = 1, mask_size_l(mid,2)
    13741381                DO  k = 1, mask_size_l(mid,3)
    1375                    local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),  &
    1376                                                            mask_j(mid,j),  &
    1377                                                            mask_i(mid,i)), &
    1378                                             REAL( fill_value, KIND = wp ), &
    1379                                             BTEST( wall_flags_total_0(     &
    1380                                                            mask_k(mid,k),  &
    1381                                                            mask_j(mid,j),  &
    1382                                                            mask_i(mid,i)), &
     1382                   local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),                          &
     1383                                                           mask_j(mid,j),                          &
     1384                                                           mask_i(mid,i)),                         &
     1385                                            REAL( fill_value, KIND = wp ),                         &
     1386                                            BTEST( wall_flags_total_0(mask_k(mid,k),               &
     1387                                                                      mask_j(mid,j),               &
     1388                                                                      mask_i(mid,i)),              &
    13831389                                                   flag_nr ) )
    13841390                ENDDO
    … …  
    13951401                im = mask_i(mid,i)
    13961402                jm = mask_j(mid,j)
    1397                 ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
    1398                               DIM = 1 ) - 1
     1403                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1
    13991404                DO  k = 1, mask_size_l(mid,3)
    14001405                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
    … …  
    14151420 END SUBROUTINE doq_output_mask
    14161421
    1417 !------------------------------------------------------------------------------!
     1422!--------------------------------------------------------------------------------------------------!
    14181423! Description:
    14191424! ------------
    14201425!> Allocate required arrays
    1421 !------------------------------------------------------------------------------!
     1426!--------------------------------------------------------------------------------------------------!
    14221427 SUBROUTINE doq_init
    14231428
    … …  
    14251430
    14261431    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
     1432
    14271433
    14281434!
    … …  
    15641570
    15651571!
    1566 !-- Save timestep number to check in time_integration if doq_calculate
    1567 !-- has been called already, since the CALL occurs at two locations, but the calculations need to be
    1568 !-- done only once per timestep.
     1572!-- Save timestep number to check in time_integration if doq_calculate has been called already,
     1573!-- since the CALL occurs at two locations, but the calculations need to be done only once per
     1574!-- timestep.
    15691575    timestep_number_at_prev_calc = current_timestep_number
    15701576
    … …  
    15801586                DO  j = nys, nyn
    15811587                   DO  k = nzb+1, nzt
    1582                       ti(k,j,i) = 0.25_wp * SQRT(                              &
    1583                         (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
    1584                               - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
    1585                           - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
    1586                               - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
    1587                       + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
    1588                               - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
    1589                           - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
    1590                               - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
    1591                       + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
    1592                               - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
    1593                           - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
    1594                               - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
    1595                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
     1588                      ti(k,j,i) = 0.25_wp * SQRT(                                                  &
     1589                                       (   (   w(k,j+1,i) + w(k-1,j+1,i)                           &
     1590                                             - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                   &
     1591                                         - (   v(k+1,j,i) + v(k+1,j+1,i)                           &
     1592                                             - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2          &
     1593                                     + (   (   u(k+1,j,i) + u(k+1,j,i+1)                           &
     1594                                             - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)               &
     1595                                         - (   w(k,j,i+1) + w(k-1,j,i+1)                           &
     1596                                             - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2          &
     1597                                     + (   (   v(k,j,i+1) + v(k,j+1,i+1)                           &
     1598                                             - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                   &
     1599                                         - (   u(k,j+1,i) + u(k,j+1,i+1)                           &
     1600                                             - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )       &
     1601                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    15961602                   ENDDO
    15971603                ENDDO
    … …  
    16031609                DO  j = nys, nyn
    16041610                   DO  k = nzb+1, nzt
    1605                       uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
    1606                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1607                                 BTEST( wall_flags_total_0(k,j,i), 1) )
     1611                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                                              &
     1612                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    16081613                   ENDDO
    16091614                ENDDO
    … …  
    16151620                DO  j = nys, nyn
    16161621                   DO  k = nzb+1, nzt
    1617                       vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
    1618                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1619                                 BTEST( wall_flags_total_0(k,j,i), 2) )
     1622                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                                              &
     1623                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    16201624                   ENDDO
    16211625                ENDDO
    … …  
    16271631                DO  j = nys, nyn
    16281632                   DO  k = nzb+1, nzt-1
    1629                       ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
    1630                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1631                                 BTEST( wall_flags_total_0(k,j,i), 3) )
     1633                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                                              &
     1634                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    16321635                   ENDDO
    16331636                ENDDO
    … …  
    16391642                DO  j = nys, nyn
    16401643                   DO  k = nzb+1, nzt-1
    1641                       wu(k,j,i) = w(k,j,i)                                     &
    1642                        * 0.25_wp * ( u(k,j,i)   + u(k,j,i+1)                   &
    1643                                    + u(k+1,j,i) + u(k+1,j,i+1) )               &
    1644                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1645                                 BTEST( wall_flags_total_0(k,j,i), 0) )
     1644                      wu(k,j,i) = w(k,j,i)                                                         &
     1645                                  * 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )&
     1646                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    16461647                   ENDDO
    16471648                ENDDO
    … …  
    16531654                DO  j = nys, nyn
    16541655                   DO  k = nzb+1, nzt-1
    1655                       wv(k,j,i) = w(k,j,i)                                     &
    1656                        * 0.25_wp * ( v(k,j,i)   + v(k,j+1,i)                   &
    1657                                    + v(k+1,j,i) + v(k+1,j+1,i) )               &
    1658                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1659                                 BTEST( wall_flags_total_0(k,j,i), 0) )
     1656                      wv(k,j,i) = w(k,j,i)                                                         &
     1657                                  * 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )&
     1658                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    16601659                   ENDDO
    16611660                ENDDO
    … …  
    16671666                DO  j = nys, nyn
    16681667                   DO  k = nzb+1, nzt-1
    1669                       wtheta(k,j,i) = w(k,j,i)                                 &
    1670                        *  0.5_wp  * ( pt(k,j,i) + pt(k+1,j,i) )                &
    1671                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1672                                 BTEST( wall_flags_total_0(k,j,i), 3) )
     1668                      wtheta(k,j,i) = w(k,j,i)                                                     &
     1669                                     *  0.5_wp  * ( pt(k,j,i) + pt(k+1,j,i) )                      &
     1670                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ))
    16731671                   ENDDO
    16741672                ENDDO
    … …  
    16801678                DO  j = nys, nyn
    16811679                   DO  k = nzb+1, nzt-1
    1682                       wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )&
    1683                        * MERGE( 1.0_wp, 0.0_wp,                                &
    1684                                 BTEST( wall_flags_total_0(k,j,i), 3) )
     1680                      wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )                    &
     1681                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    16851682                   ENDDO
    16861683                ENDDO
    … …  
    16901687          CASE ( 'theta_2m*' )
    16911688!
    1692 !--          2-m potential temperature is caluclated from surface arrays. In
    1693 !--          case the 2m level is below the first grid point, MOST is applied,
    1694 !--          else, linear interpolation between two vertical grid levels is
    1695 !--          applied. To access all surfaces, iterate over all horizontally-
     1689!--          2-m potential temperature is caluclated from surface arrays. In case the 2m level is
     1690!--          below the first grid point, MOST is applied, else, linear interpolation between two
     1691!--          vertical grid levels is applied. To access all surfaces, iterate over all horizontally-
    16961692!--          upward facing surface types.
    16971693             surf => surf_def_h(0)
    … …  
    17051701          CASE ( 'wspeed_10m*' )
    17061702!
    1707 !--          10-m wind speed is caluclated from surface arrays. In
    1708 !--          case the 10m level is below the first grid point, MOST is applied,
    1709 !--          else, linear interpolation between two vertical grid levels is
    1710 !--          applied. To access all surfaces, iterate over all horizontally-
    1711 !--          upward facing surface types.
     1703!--          10-m wind speed is caluclated from surface arrays. In case the 10m level is below the
     1704!--          first grid point, MOST is applied, else, linear interpolation between two vertical grid
     1705!--          levels is applied. To access all surfaces, iterate over all horizontally-upward facing
     1706!--          surface types.
    17121707             surf => surf_def_h(0)
    17131708             CALL calc_wind_10m
    … …  
    17171712             CALL calc_wind_10m
    17181713!
    1719 !--       horizontal wind speed
     1714!--       Horizontal wind speed
    17201715          CASE ( 'wspeed' )
    17211716             DO  i = nxl, nxr
    17221717                DO  j = nys, nyn
    17231718                   DO  k = nzb, nzt+1
    1724                       wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2             &
    1725                                           + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )           &
    1726                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
     1719                      wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2              &
     1720                                          + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )            &
     1721                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ))
    17271722                   ENDDO
    17281723                ENDDO
    … …  
    17301725
    17311726!
    1732 !--       horizontal wind direction
     1727!--       Horizontal wind direction
    17331728          CASE ( 'wdir' )
    17341729             DO  i = nxl, nxr
    … …  
    17381733                      v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
    17391734
    1740                       wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &
    1741                                   / pi * 180.0_wp + 180.0_wp
     1735                      wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) )                      &
     1736                                    / pi * 180.0_wp + 180.0_wp
    17421737                   ENDDO
    17431738                ENDDO
    … …  
    17521747
    17531748!
    1754 !-- The following block contains subroutines to calculate diagnostic
    1755 !-- surface quantities.
     1749!-- The following block contains subroutines to calculate diagnostic surface quantities.
    17561750    CONTAINS
    1757 !------------------------------------------------------------------------------!
     1751!--------------------------------------------------------------------------------------------------!
    17581752! Description:
    17591753! ------------
    17601754!> Calculation of 2-m potential temperature.
    1761 !------------------------------------------------------------------------------!
     1755!--------------------------------------------------------------------------------------------------!
    17621756       SUBROUTINE calc_pt_2m
    17631757
    1764           USE surface_layer_fluxes_mod,                                        &
     1758          USE surface_layer_fluxes_mod,                                                            &
    17651759              ONLY:  psi_h
    17661760
    … …  
    17761770             k = surf%k(m)
    17771771!
    1778 !--          If 2-m level is below the first grid level, MOST is
    1779 !--          used for calculation of 2-m temperature.
     1772!--          If 2-m level is below the first grid level, MOST is used for calculation of
     1773!--          2-m temperature.
    17801774             IF ( surf%z_mo(m) > 2.0_wp )  THEN
    1781                 pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa           &
    1782                                 * ( LOG( 2.0_wp /  surf%z0h(m) )               &
    1783                                   - psi_h( 2.0_wp / surf%ol(m) )               &
    1784                                   + psi_h( surf%z0h(m) / surf%ol(m) ) )
    1785 !
    1786 !--          If 2-m level is above the first grid level, 2-m temperature
    1787 !--          is linearly interpolated between the two nearest vertical grid
    1788 !--          levels. Note, since 2-m temperature is only computed for
    1789 !--          horizontal upward-facing surfaces, only a vertical
    1790 !--          interpolation is necessary.
     1775                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa                               &
     1776                                * ( LOG( 2.0_wp / surf%z0h(m) )                                    &
     1777                                    - psi_h( 2.0_wp      / surf%ol(m) )                            &
     1778                                    + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1779!
     1780!--          If 2-m level is above the first grid level, 2-m temperature is linearly interpolated
     1781!--          between the two nearest vertical grid levels. Note, since 2-m temperature is only
     1782!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
     1783!--          necessary.
    17911784             ELSE
    17921785!
    … …  
    17981791!
    17991792!--             kk defines the index of the first grid level >= 2m.
    1800                 pt_2m(j,i) = pt(kk-1,j,i) +                                    &
    1801                               ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
    1802                               ( pt(kk,j,i)       - pt(kk-1,j,i) ) /            &
     1793                pt_2m(j,i) = pt(kk-1,j,i) +                                                        &
     1794                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *                                &
     1795                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /                                &
    18031796                              ( zu(kk)           - zu(kk-1)     )
    18041797             ENDIF
    … …  
    18081801       END SUBROUTINE calc_pt_2m
    18091802
    1810 !------------------------------------------------------------------------------!
     1803!--------------------------------------------------------------------------------------------------!
    18111804! Description:
    18121805! ------------
    18131806!> Calculation of 10-m wind speed.
    1814 !------------------------------------------------------------------------------!
     1807!--------------------------------------------------------------------------------------------------!
    18151808       SUBROUTINE calc_wind_10m
    18161809
    1817           USE surface_layer_fluxes_mod,                                        &
     1810          USE surface_layer_fluxes_mod,                                                            &
    18181811              ONLY:  psi_m
    18191812
    … …  
    18251818          REAL(wp) ::  uv_l !< wind speed at lower grid point
    18261819          REAL(wp) ::  uv_u !< wind speed at upper grid point
     1820
    18271821
    18281822          DO  m = 1, surf%ns
    … …  
    18321826             k = surf%k(m)
    18331827!
    1834 !--          If 10-m level is below the first grid level, MOST is
    1835 !--          used for calculation of 10-m temperature.
     1828!--          If 10-m level is below the first grid level, MOST is used for calculation of 10-m
     1829!--          temperature.
    18361830             IF ( surf%z_mo(m) > 10.0_wp )  THEN
    1837                 uv_10m(j,i) = surf%us(m) / kappa                               &
    1838                           * ( LOG( 10.0_wp /  surf%z0(m) )                     &
    1839                               - psi_m( 10.0_wp    / surf%ol(m) )               &
    1840                               + psi_m( surf%z0(m) / surf%ol(m) ) )
    1841 !
    1842 !--          If 10-m level is above the first grid level, 10-m wind speed
    1843 !--          is linearly interpolated between the two nearest vertical grid
    1844 !--          levels. Note, since 10-m temperature is only computed for
    1845 !--          horizontal upward-facing surfaces, only a vertical
    1846 !--          interpolation is necessary.
     1831                uv_10m(j,i) = surf%us(m) / kappa                                                   &
     1832                              * ( LOG( 10.0_wp /  surf%z0(m) )                                     &
     1833                                  - psi_m( 10.0_wp    / surf%ol(m) )                               &
     1834                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
     1835!
     1836!--          If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated
     1837!--          between the two nearest vertical grid levels. Note, since 10-m temperature is only
     1838!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
     1839!--          necessary.
    18471840             ELSE
    18481841!
    … …  
    18541847!
    18551848!--             kk defines the index of the first grid level >= 10m.
    1856                 uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2   &
     1849                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2                       &
    18571850                           + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
    18581851
    1859                 uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2   &
     1852                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2                       &
    18601853                           + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
    18611854
    1862                 uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *        &
    1863                                      ( uv_u              - uv_l     ) /        &
     1855                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *                            &
     1856                                     ( uv_u              - uv_l     ) /                            &
    18641857                                     ( zu(kk)            - zu(kk-1) )
    18651858
    … …  
    18731866
    18741867
    1875 !------------------------------------------------------------------------------!
     1868!--------------------------------------------------------------------------------------------------!
    18761869! Description:
    18771870! ------------
    1878 !> Preparation of the diagnostic output, counting of the module-specific
    1879 !> output quantities and gathering of the output names.
    1880 !------------------------------------------------------------------------------!
     1871!> Preparation of the diagnostic output, counting of the module-specific output quantities and
     1872!> gathering of the output names.
     1873!--------------------------------------------------------------------------------------------------!
    18811874 SUBROUTINE doq_prepare
    18821875
    1883     USE control_parameters,                                                    &
     1876    USE control_parameters,                                                                        &
    18841877        ONLY:  do2d, do3d, domask, masks
    18851878
    18861879    IMPLICIT NONE
    18871880
    1888     CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !<
    1889                                                           !< label array for 2d output quantities
     1881    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !< label array for 2d output quantities
    18901882
    18911883    INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
    … …  
    18951887    INTEGER(iwp) ::  mid          !< masked output running index
    18961888
     1889
    18971890    prepared_diagnostic_output_quantities = .FALSE.
    18981891
    … …  
    19071900!
    19081901!--    Gather 2d output quantity names.
    1909 !--    Check for double occurrence of output quantity, e.g. by _xy,
    1910 !--    _yz, _xz.
     1902!--    Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz.
    19111903       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
    19121904          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
    … …  
    19301922       ivar = 1
    19311923!
    1932 !--    Gather masked output quantity names. Also check for double output
    1933 !--    e.g. by different masks.
     1924!--    Gather masked output quantity names. Also check for double output e.g. by different masks.
    19341925       DO  mid = 1, masks
    19351926          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
    … …  
    19501941 END SUBROUTINE doq_prepare
    19511942
    1952 !------------------------------------------------------------------------------!
     1943!--------------------------------------------------------------------------------------------------!
    19531944! Description:
    19541945! ------------
    19551946!> Subroutine reads local (subdomain) restart data
    1956 !------------------------------------------------------------------------------!
    1957  SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,                             &
    1958                                nxr_on_file, nynf, nync, nyn_on_file, nysf,                         &
    1959                                nysc, nys_on_file, tmp_2d, tmp_3d, found )
     1947!--------------------------------------------------------------------------------------------------!
     1948 SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,    &
     1949                               nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
    19601950
    19611951    USE control_parameters
    … …  
    20021992          pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                                     &
    20031993                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    2004    
     1994
    20051995       CASE ( 'ti_av' )
    20061996          IF ( .NOT. ALLOCATED( ti_av ) )  THEN
    … …  
    21082098 END SUBROUTINE doq_rrd_local_ftn
    21092099
    2110 !------------------------------------------------------------------------------!
     2100!--------------------------------------------------------------------------------------------------!
    21112101! Description:
    21122102! ------------
    21132103!> Read module-specific local restart data arrays (MPI-IO).
    2114 !------------------------------------------------------------------------------!
     2104!--------------------------------------------------------------------------------------------------!
    21152105 SUBROUTINE doq_rrd_local_mpi
    21162106
    … …  
    21982188 END SUBROUTINE doq_rrd_local_mpi
    21992189
    2200 !------------------------------------------------------------------------------!
     2190!--------------------------------------------------------------------------------------------------!
    22012191! Description:
    22022192! ------------
    22032193!> Subroutine writes local (subdomain) restart data
    2204 !------------------------------------------------------------------------------!
     2194!--------------------------------------------------------------------------------------------------!
    22052195 SUBROUTINE doq_wrd_local
    22062196
    2207 
    22082197    IMPLICIT NONE
     2198
    22092199
    22102200    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN

  • palm/trunk/SOURCE/diffusion_s.f90

    r4360 r4583  
    11!> @file diffusion_s.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
    27 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    28 ! topography information used in wall_flags_static_0
    29 !
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
     29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     30! information used in wall_flags_static_0
     31!
    3032! 4329 2019-12-10 15:46:36Z motisi
    3133! Renamed wall_flags_0 to wall_flags_total_0
    32 ! 
     34!
    3335! 4182 2019-08-22 15:20:23Z scharf
    3436! Corrected "Former revisions" section
    35 ! 
     37!
    3638! 3927 2019-04-23 13:24:29Z raasch
    3739! pointer attribute removed from scalar 3d-array for performance reasons
    38 ! 
     40!
    3941! 3761 2019-02-25 15:31:42Z raasch
    4042! unused variables removed
    41 ! 
     43!
    4244! 3655 2019-01-07 16:51:22Z knoop
    4345! nopointer option removed
    … …  
    5052! ------------
    5153!> Diffusion term of scalar quantities (temperature and water content)
    52 !------------------------------------------------------------------------------!
     54!--------------------------------------------------------------------------------------------------!
    5355 MODULE diffusion_s_mod
    54  
     56
    5557
    5658    PRIVATE
    … …  
    6567
    6668
    67 !------------------------------------------------------------------------------!
     69!--------------------------------------------------------------------------------------------------!
    6870! Description:
    6971! ------------
    7072!> Call for all grid points
    71 !------------------------------------------------------------------------------!
    72     SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,          &
    73                                s_flux_t,                                       &
    74                                s_flux_lsm_h_up,    s_flux_usm_h_up,            &
    75                                s_flux_def_v_north, s_flux_def_v_south,         &
    76                                s_flux_def_v_east,  s_flux_def_v_west,          &
    77                                s_flux_lsm_v_north, s_flux_lsm_v_south,         &
    78                                s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
    79                                s_flux_usm_v_north, s_flux_usm_v_south,         &
     73!--------------------------------------------------------------------------------------------------!
     74    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,                              &
     75                               s_flux_t,                                                           &
     76                               s_flux_lsm_h_up,    s_flux_usm_h_up,                                &
     77                               s_flux_def_v_north, s_flux_def_v_south,                             &
     78                               s_flux_def_v_east,  s_flux_def_v_west,                              &
     79                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
     80                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
     81                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
    8082                               s_flux_usm_v_east,  s_flux_usm_v_west )
    8183
    82        USE arrays_3d,                                                          &
     84       USE arrays_3d,                                                                              &
    8385           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    84        
    85        USE control_parameters,                                                 &
     86
     87       USE control_parameters,                                                                     &
    8688           ONLY: use_surface_fluxes, use_top_fluxes
    87        
    88        USE grid_variables,                                                     &
     89
     90       USE grid_variables,                                                                         &
    8991           ONLY:  ddx, ddx2, ddy, ddy2
    90        
    91        USE indices,                                                            &
    92            ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,        &
    93                   wall_flags_total_0
    94        
     92
     93       USE indices,                                                                                &
     94           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
     95
    9596       USE kinds
    9697
    97        USE surface_mod,                                                        &
    98            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    99                    surf_usm_v
     98       USE surface_mod,                                                                            &
     99           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    100100
    101101       IMPLICIT NONE
    … …  
    109109
    110110       REAL(wp) ::  flag              !< flag to mask topography grid points
    111        REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
    112        REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
     111       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
     112       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
    113113       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
    114        REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     114       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     115       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
    115116       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
    116        REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
    117 
     117
     118       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
     119       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
     120       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
    118121       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
    119122       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
    120        REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
    121123       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
    122        REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
    123        REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
    124124       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
     125       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
    125126       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
    126127       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
    127        REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
    128128       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
    129129       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
     130       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
    130131       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
    131132       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
    132        REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
    133133       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
    134134       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
    … …  
    164164!
    165165!--             Predetermine flag to mask topography and wall-bounded grid points
    166                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 
    167 !
    168 !--             Predetermine flag to mask wall-bounded grid points, equivalent to
    169 !--             former s_outer array
     166                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     167!
     168!--             Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer
     169!--             array
    170170                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
    171171                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
    … …  
    173173                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
    174174
    175                 tend(k,j,i) = tend(k,j,i)                                      &
    176                                           + 0.5_wp * (                         &
    177                         mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )               &
    178                                    * ( s(k,j,i+1) - s(k,j,i)   )               &
    179                       - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )               &
    180                                    * ( s(k,j,i)   - s(k,j,i-1) )               &
    181                                                      ) * ddx2 * flag           &
    182                                           + 0.5_wp * (                         &
    183                         mask_north * ( kh(k,j,i) + kh(k,j+1,i) )               &
    184                                    * ( s(k,j+1,i) - s(k,j,i)   )               &
    185                       - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )               &
    186                                    * ( s(k,j,i)   - s(k,j-1,i) )               &
    187                                                      ) * ddy2 * flag
    188              ENDDO
    189 
    190 !
    191 !--          Apply prescribed horizontal wall heatflux where necessary. First,
    192 !--          determine start and end index for respective (j,i)-index. Please
    193 !--          note, in the flat case following loop will not be entered, as
    194 !--          surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
    195 !--          so far.
    196 !--          First, for default-type surfaces
     175                tend(k,j,i) = tend(k,j,i)                                                          &
     176                                                  + 0.5_wp * (                                     &
     177                                mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                           &
     178                                           * ( s(k,j,i+1) - s(k,j,i)   )                           &
     179                              - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                           &
     180                                           * ( s(k,j,i)   - s(k,j,i-1) )                           &
     181                                                             ) * ddx2 * flag                       &
     182                                                  + 0.5_wp * (                                     &
     183                                mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                           &
     184                                           * ( s(k,j+1,i) - s(k,j,i)   )                           &
     185                              - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                           &
     186                                           * ( s(k,j,i)   - s(k,j-1,i) )                           &
     187                                                             ) * ddy2 * flag
     188             ENDDO
     189
     190!
     191!--          Apply prescribed horizontal wall heatflux where necessary. First, determine start and
     192!--          end index for respective (j,i)-index. Please note, in the flat case following loop will
     193!--          not be entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural
     194!--          surfaces so far.
     195!--          First, for default-type surfaces.
    197196!--          North-facing vertical default-type surfaces
    198197             surf_s = surf_def_v(0)%start_index(j,i)
    … …  
    294293
    295294!
    296 !--          Compute vertical diffusion. In case that surface fluxes have been
    297 !--          prescribed or computed at bottom and/or top, index k starts/ends at
    298 !--          nzb+2 or nzt-1, respectively. Model top is also mask if top flux
    299 !--          is given.
     295!--          Compute vertical diffusion. In case that surface fluxes have been prescribed or
     296!--          computed at bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively.
     297!--          Model top is also mask if top flux is given.
    300298             DO  k = nzb+1, nzt
    301299!
    302 !--             Determine flags to mask topography below and above. Flag 0 is
    303 !--             used to mask topography in general, and flag 8 implies
    304 !--             information about use_surface_fluxes. Flag 9 is used to control
    305 !--             flux at model top.
    306                 mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
    307                                  BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    308                 mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
    309                                  BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
    310                               MERGE( 1.0_wp, 0.0_wp,                           &
    311                                  BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    312                 flag        = MERGE( 1.0_wp, 0.0_wp,                           &
    313                                  BTEST( wall_flags_total_0(k,j,i), 0 ) )
    314 
    315                 tend(k,j,i) = tend(k,j,i)                                      &
    316                                        + 0.5_wp * (                            &
    317                                       ( kh(k,j,i) + kh(k+1,j,i) ) *            &
    318                                           ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    319                                                             * rho_air_zw(k)    &
    320                                                             * mask_top         &
    321                                     - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
    322                                           ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    323                                                             * rho_air_zw(k-1)  &
    324                                                             * mask_bottom      &
    325                                                   ) * ddzw(k) * drho_air(k)    &
     300!--             Determine flags to mask topography below and above. Flag 0 is used to mask
     301!--             topography in general, and flag 8 implies information about use_surface_fluxes.
     302!--             Flag 9 is used to control flux at model top.
     303                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     304                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
     305                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     306                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     307
     308                tend(k,j,i) = tend(k,j,i)                                                          &
     309                                       + 0.5_wp * (                                                &
     310                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
     311                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
     312                                                            * rho_air_zw(k)                        &
     313                                                            * mask_top                             &
     314                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
     315                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
     316                                                            * rho_air_zw(k-1)                      &
     317                                                            * mask_bottom                          &
     318                                                  ) * ddzw(k) * drho_air(k)                        &
    326319                                                              * flag
    327320             ENDDO
    … …  
    331324             IF ( use_surface_fluxes )  THEN
    332325!
    333 !--             Default-type surfaces, upward-facing               
     326!--             Default-type surfaces, upward-facing
    334327                surf_s = surf_def_h(0)%start_index(j,i)
    335328                surf_e = surf_def_h(0)%end_index(j,i)
    … …  
    337330
    338331                   k   = surf_def_h(0)%k(m)
    339                    tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)              &
    340                                        * ddzw(k) * drho_air(k)
     332                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
    341333
    342334                ENDDO
    343335!
    344 !--             Default-type surfaces, downward-facing               
     336!--             Default-type surfaces, downward-facing
    345337                surf_s = surf_def_h(1)%start_index(j,i)
    346338                surf_e = surf_def_h(1)%end_index(j,i)
    … …  
    348340
    349341                   k   = surf_def_h(1)%k(m)
    350                    tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)            &
    351                                        * ddzw(k) * drho_air(k)
     342                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
    352343
    353344                ENDDO
    354345!
    355 !--             Natural-type surfaces, upward-facing 
     346!--             Natural-type surfaces, upward-facing
    356347                surf_s = surf_lsm_h%start_index(j,i)
    357348                surf_e = surf_lsm_h%end_index(j,i)
    … …  
    359350
    360351                   k   = surf_lsm_h%k(m)
    361                    tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)              &
    362                                        * ddzw(k) * drho_air(k)
     352                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
    363353
    364354                ENDDO
    365355!
    366 !--             Urban-type surfaces, upward-facing     
     356!--             Urban-type surfaces, upward-facing
    367357                surf_s = surf_usm_h%start_index(j,i)
    368358                surf_e = surf_usm_h%end_index(j,i)
    … …  
    370360
    371361                   k   = surf_usm_h%k(m)
    372                    tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)              &
    373                                        * ddzw(k) * drho_air(k)
     362                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
    374363
    375364                ENDDO
    … …  
    384373
    385374                   k   = surf_def_h(2)%k(m)
    386                    tend(k,j,i) = tend(k,j,i)                                   &
    387                            + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
     375                   tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
    388376                ENDDO
    389377             ENDIF
    … …  
    394382    END SUBROUTINE diffusion_s
    395383
    396 !------------------------------------------------------------------------------!
     384!--------------------------------------------------------------------------------------------------!
    397385! Description:
    398386! ------------
    399387!> Call for grid point i,j
    400 !------------------------------------------------------------------------------!
    401     SUBROUTINE diffusion_s_ij( i, j, s,                                        &
    402                                s_flux_def_h_up,    s_flux_def_h_down,          &
    403                                s_flux_t,                                       &
    404                                s_flux_lsm_h_up,    s_flux_usm_h_up,            &
    405                                s_flux_def_v_north, s_flux_def_v_south,         &
    406                                s_flux_def_v_east,  s_flux_def_v_west,          &
    407                                s_flux_lsm_v_north, s_flux_lsm_v_south,         &
    408                                s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
    409                                s_flux_usm_v_north, s_flux_usm_v_south,         &
    410                                s_flux_usm_v_east,  s_flux_usm_v_west )       
    411 
    412        USE arrays_3d,                                                          &
     388!--------------------------------------------------------------------------------------------------!
     389    SUBROUTINE diffusion_s_ij( i, j, s,                                                            &
     390                               s_flux_def_h_up,    s_flux_def_h_down,                              &
     391                               s_flux_t,                                                           &
     392                               s_flux_lsm_h_up,    s_flux_usm_h_up,                                &
     393                               s_flux_def_v_north, s_flux_def_v_south,                             &
     394                               s_flux_def_v_east,  s_flux_def_v_west,                              &
     395                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
     396                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
     397                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
     398                               s_flux_usm_v_east,  s_flux_usm_v_west )
     399
     400       USE arrays_3d,                                                                              &
    413401           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
    414            
    415        USE control_parameters,                                                 &
     402
     403       USE control_parameters,                                                                     &
    416404           ONLY: use_surface_fluxes, use_top_fluxes
    417        
    418        USE grid_variables,                                                     &
     405
     406       USE grid_variables,                                                                         &
    419407           ONLY:  ddx, ddx2, ddy, ddy2
    420        
    421        USE indices,                                                            &
     408
     409       USE indices,                                                                                &
    422410           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_total_0
    423        
     411
    424412       USE kinds
    425413
    426        USE surface_mod,                                                        &
    427            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    428                    surf_usm_v
     414       USE surface_mod,                                                                            &
     415           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    429416
    430417       IMPLICIT NONE
    … …  
    438425
    439426       REAL(wp) ::  flag              !< flag to mask topography grid points
    440        REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
    441        REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point 
     427       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
     428       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
    442429       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
    443        REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     430       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
     431       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
    444432       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
    445        REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
    446 
     433
     434       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
     435       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
     436       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
    447437       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
    448438       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
    449        REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
    450439       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
    451        REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
    452        REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
    453440       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
     441       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
    454442       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
    455443       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
    456        REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
    457444       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
    458445       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
     446       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
    459447       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
    460448       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
    461        REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
    462449       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
    463450       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
    … …  
    465452       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< treated scalar
    466453
     454
    467455!
    468456!--    Compute horizontal diffusion
    … …  
    470458!
    471459!--       Predetermine flag to mask topography and wall-bounded grid points
    472           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    473 !
    474 !--       Predetermine flag to mask wall-bounded grid points, equivalent to
    475 !--       former s_outer array
     460          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     461!
     462!--       Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer array
    476463          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
    477464          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
    … …  
    479466          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
    480467!
    481 !--       Finally, determine flag to mask both topography itself as well
    482 !--       as wall-bounded grid points, which will be treated further below
    483 
    484           tend(k,j,i) = tend(k,j,i)                                            &
    485                                           + 0.5_wp * (                         &
    486                             mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )           &
    487                                        * ( s(k,j,i+1) - s(k,j,i)   )           &
    488                           - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )           &
    489                                        * ( s(k,j,i)   - s(k,j,i-1) )           &
    490                                                      ) * ddx2 * flag           &
    491                                           + 0.5_wp * (                         &
    492                             mask_north * ( kh(k,j,i) + kh(k,j+1,i) )           &
    493                                        * ( s(k,j+1,i) - s(k,j,i)   )           &
    494                           - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )           &
    495                                        * ( s(k,j,i)  - s(k,j-1,i)  )           &
     468!--       Finally, determine flag to mask both topography itself as well as wall-bounded grid
     469!--       points, which will be treated further below
     470
     471          tend(k,j,i) = tend(k,j,i)                                                                &
     472                                          + 0.5_wp * (                                             &
     473                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                               &
     474                                       * ( s(k,j,i+1) - s(k,j,i)   )                               &
     475                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                               &
     476                                       * ( s(k,j,i)   - s(k,j,i-1) )                               &
     477                                                     ) * ddx2 * flag                               &
     478                                          + 0.5_wp * (                                             &
     479                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                               &
     480                                       * ( s(k,j+1,i) - s(k,j,i)   )                               &
     481                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                               &
     482                                       * ( s(k,j,i)  - s(k,j-1,i)  )                               &
    496483                                                     ) * ddy2 * flag
    497484       ENDDO
    498485
    499486!
    500 !--    Apply prescribed horizontal wall heatflux where necessary. First,
    501 !--    determine start and end index for respective (j,i)-index. Please
    502 !--    note, in the flat case following loops will not be entered, as
    503 !--    surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
    504 !--    so far.
    505 !--    First, for default-type surfaces
     487!--    Apply prescribed horizontal wall heatflux where necessary. First, determine start and end
     488!--    index for respective (j,i)-index. Please note, in the flat case following loops will not be
     489!--    entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces so far.
     490!--    First, for default-type surfaces.
    506491!--    North-facing vertical default-type surfaces
    507492       surf_s = surf_def_v(0)%start_index(j,i)
    … …  
    604589
    605590!
    606 !--    Compute vertical diffusion. In case that surface fluxes have been
    607 !--    prescribed or computed at bottom and/or top, index k starts/ends at
    608 !--    nzb+2 or nzt-1, respectively. Model top is also mask if top flux
    609 !--    is given.
     591!--    Compute vertical diffusion. In case that surface fluxes have been prescribed or computed at
     592!--    bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively. Model top is also
     593!--    mask if top flux is given.
    610594       DO  k = nzb+1, nzt
    611595!
    612 !--       Determine flags to mask topography below and above. Flag 0 is
    613 !--       used to mask topography in general, and flag 8 implies
    614 !--       information about use_surface_fluxes. Flag 9 is used to control
    615 !--       flux at model top.   
    616           mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
    617                                BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    618           mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
    619                                BTEST( wall_flags_total_0(k+1,j,i), 8 ) )  *    &
    620                         MERGE( 1.0_wp, 0.0_wp,                                 &
    621                                BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    622           flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
    623                                BTEST( wall_flags_total_0(k,j,i), 0 ) )
    624 
    625           tend(k,j,i) = tend(k,j,i)                                            &
    626                                        + 0.5_wp * (                            &
    627                                       ( kh(k,j,i) + kh(k+1,j,i) ) *            &
    628                                           ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
    629                                                             * rho_air_zw(k)    &
    630                                                             * mask_top         &
    631                                     - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
    632                                           ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
    633                                                             * rho_air_zw(k-1)  &
    634                                                             * mask_bottom      &
    635                                                   ) * ddzw(k) * drho_air(k)    &
     596!--       Determine flags to mask topography below and above. Flag 0 is used to mask topography in
     597!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
     598!--       control flux at model top.
     599          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     600          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) )  *        &
     601                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     602          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     603
     604          tend(k,j,i) = tend(k,j,i)                                                                &
     605                                       + 0.5_wp * (                                                &
     606                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
     607                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
     608                                                            * rho_air_zw(k)                        &
     609                                                            * mask_top                             &
     610                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
     611                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
     612                                                            * rho_air_zw(k-1)                      &
     613                                                            * mask_bottom                          &
     614                                                  ) * ddzw(k) * drho_air(k)                        &
    636615                                                              * flag
    637616       ENDDO
    … …  
    649628             k   = surf_def_h(0)%k(m)
    650629
    651              tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)                    &
    652                                        * ddzw(k) * drho_air(k)
     630             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
    653631          ENDDO
    654632!
    … …  
    660638             k   = surf_def_h(1)%k(m)
    661639
    662              tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)                  &
    663                                        * ddzw(k) * drho_air(k)
     640             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
    664641          ENDDO
    665642!
    … …  
    670647             k   = surf_lsm_h%k(m)
    671648
    672              tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)                    &
    673                                        * ddzw(k) * drho_air(k)
     649             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
    674650          ENDDO
    675651!
    … …  
    680656             k   = surf_usm_h%k(m)
    681657
    682              tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)                    &
    683                                        * ddzw(k) * drho_air(k)
     658             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
    684659          ENDDO
    685660       ENDIF
    … …  
    692667
    693668             k   = surf_def_h(2)%k(m)
    694              tend(k,j,i) = tend(k,j,i)                                         &
    695                            + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
     669             tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
    696670          ENDDO
    697671       ENDIF

  • palm/trunk/SOURCE/diffusion_u.f90

    r4360 r4583  
    11!> @file diffusion_u.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
    27 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    28 ! topography information used in wall_flags_static_0
    29 !
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
     29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     30! information used in wall_flags_static_0
     31!
    3032! 4329 2019-12-10 15:46:36Z motisi
    3133! Renamed wall_flags_0 to wall_flags_static_0
    32 ! 
     34!
    3335! 4182 2019-08-22 15:20:23Z scharf
    3436! Corrected "Former revisions" section
    35 ! 
     37!
    3638! 3655 2019-01-07 16:51:22Z knoop
    3739! OpenACC port for SPEC
    … …  
    4446! ------------
    4547!> Diffusion term of the u-component
    46 !> @todo additional damping (needed for non-cyclic bc) causes bad vectorization
    47 !>       and slows down the speed on NEC about 5-10%
    48 !------------------------------------------------------------------------------!
     48!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the
     49!        speed on NEC about 5-10%
     50!--------------------------------------------------------------------------------------------------!
    4951 MODULE diffusion_u_mod
    50  
     52
    5153
    5254    PRIVATE
    … …  
    6163
    6264
    63 !------------------------------------------------------------------------------!
     65!--------------------------------------------------------------------------------------------------!
    6466! Description:
    6567! ------------
    6668!> Call for all grid points
    67 !------------------------------------------------------------------------------!
     69!--------------------------------------------------------------------------------------------------!
    6870    SUBROUTINE diffusion_u
    6971
    70        USE arrays_3d,                                                          &
    71            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    72        
    73        USE control_parameters,                                                 &
    74            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    75                   use_top_fluxes
    76        
    77        USE grid_variables,                                                     &
     72       USE arrays_3d,                                                                              &
     73           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
     74
     75       USE control_parameters,                                                                     &
     76           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
     77
     78       USE grid_variables,                                                                         &
    7879           ONLY:  ddx, ddx2, ddy
    79        
    80        USE indices,                                                            &
     80
     81       USE indices,                                                                                &
    8182           ONLY:  nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
    82      
     83
    8384       USE kinds
    8485
    85        USE surface_mod,                                                        &
    86            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    87                    surf_usm_v
     86       USE surface_mod,                                                                            &
     87           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    8888
    8989       IMPLICIT NONE
    … …  
    102102       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
    103103       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
    104        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
    105        REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point 
    106        REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point 
    107        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface 
     104       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     105       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
     106       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
     107       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    108108
    109109
    … …  
    125125             DO  k = nzb+1, nzt
    126126!
    127 !--             Predetermine flag to mask topography and wall-bounded grid points. 
    128 !--             It is sufficient to masked only north- and south-facing surfaces, which
    129 !--             need special treatment for the u-component.
    130                 flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) ) 
     127!--             Predetermine flag to mask topography and wall-bounded grid points.
     128!--             It is sufficient to masked only north- and south-facing surfaces, which need special
     129!--             treatment for the u-component.
     130                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
    131131                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
    132132                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
    133133!
    134134!--             Interpolate eddy diffusivities on staggered gridpoints
    135                 kmyp = 0.25_wp *                                               &
    136                        ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
    137                 kmym = 0.25_wp *                                               &
    138                        ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    139 
    140                 tend(k,j,i) = tend(k,j,i)                                      &
    141                         + 2.0_wp * (                                           &
    142                                   km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
    143                                 - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
    144                                    ) * ddx2 * flag                             &
    145                         +          ( mask_north * (                            &
    146                             kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
    147                           + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
    148                                                   )                            &
    149                                    - mask_south * (                            &
    150                             kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    151                           + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    152                                                   )                            &
    153                                    ) * ddy  * flag                             
     135                kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
     136                kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
     137
     138                tend(k,j,i) = tend(k,j,i)                                                          &
     139                              + 2.0_wp * (                                                         &
     140                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                  &
     141                                      - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                  &
     142                                         ) * ddx2 * flag                                           &
     143                              +          ( mask_north * (                                          &
     144                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy                       &
     145                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                       &
     146                                                        )                                          &
     147                                         - mask_south * (                                          &
     148                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy                           &
     149                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx                           &
     150                                                        )                                          &
     151                                         ) * ddy  * flag
    154152             ENDDO
    155153!
    156 !--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
    157 !--          surfaces. Note, in the the flat case, loops won't be entered as
    158 !--          start_index > end_index. Furtermore, note, no vertical natural surfaces
    159 !--          so far.           
     154!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces.
     155!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
     156!--          Furtermore, note, no vertical natural surfaces so far.
    160157!--          Default-type surfaces
    161158             DO  l = 0, 1
    … …  
    164161                DO  m = surf_s, surf_e
    165162                   k           = surf_def_v(l)%k(m)
    166                    tend(k,j,i) = tend(k,j,i) +                                 &                   
    167                                     surf_def_v(l)%mom_flux_uv(m) * ddy
    168                 ENDDO   
     163                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
     164                ENDDO
    169165             ENDDO
    170166!
    … …  
    175171                DO  m = surf_s, surf_e
    176172                   k           = surf_lsm_v(l)%k(m)
    177                    tend(k,j,i) = tend(k,j,i) +                                 &                   
    178                                     surf_lsm_v(l)%mom_flux_uv(m) * ddy
    179                 ENDDO   
     173                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
     174                ENDDO
    180175             ENDDO
    181176!
    … …  
    186181                DO  m = surf_s, surf_e
    187182                   k           = surf_usm_v(l)%k(m)
    188                    tend(k,j,i) = tend(k,j,i) +                                 &                   
    189                                     surf_usm_v(l)%mom_flux_uv(m) * ddy
    190                 ENDDO   
     183                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
     184                ENDDO
    191185             ENDDO
    192186
    193187!
    194 !--          Compute vertical diffusion. In case of simulating a surface layer,
    195 !--          respective grid diffusive fluxes are masked (flag 8) within this
    196 !--          loop, and added further below, else, simple gradient approach is
    197 !--          applied. Model top is also mask if top-momentum flux is given.
     188!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
     189!--          diffusive fluxes are masked (flag 8) within this loop, and added further below, else,
     190!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
     191!--          given.
    198192             DO  k = nzb+1, nzt
    199193!
    200 !--             Determine flags to mask topography below and above. Flag 1 is
    201 !--             used to mask topography in general, and flag 8 implies
    202 !--             information about use_surface_fluxes. Flag 9 is used to control
    203 !--             momentum flux at model top. 
    204                 mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
    205                                  BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    206                 mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
    207                                  BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
    208                               MERGE( 1.0_wp, 0.0_wp,                           &
    209                                  BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    210                 flag        = MERGE( 1.0_wp, 0.0_wp,                           &
    211                                  BTEST( wall_flags_total_0(k,j,i), 1 ) )
     194!--             Determine flags to mask topography below and above. Flag 1 is used to mask
     195!--             topography in general, and flag 8 implies information about use_surface_fluxes.
     196!--             Flag 9 is used to control momentum flux at model top.
     197                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     198                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
     199                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     200                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    212201!
    213202!--             Interpolate eddy diffusivities on staggered gridpoints
    214                 kmzp = 0.25_wp *                                               &
    215                        ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
    216                 kmzm = 0.25_wp *                                               &
    217                        ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    218 
    219                 tend(k,j,i) = tend(k,j,i)                                      &
    220                         + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
    221                                    + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    222                                    ) * rho_air_zw(k)   * mask_top              &
    223                           - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
    224                                    + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
    225                                    ) * rho_air_zw(k-1) * mask_bottom           &
    226                           ) * ddzw(k) * drho_air(k) * flag
     203                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
     204                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
     205
     206                tend(k,j,i) = tend(k,j,i)                                                          &
     207                              + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                 &
     208                                         + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                       &
     209                                         ) * rho_air_zw(k)   * mask_top                            &
     210                                - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                 &
     211                                         + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                     &
     212                                         ) * rho_air_zw(k-1) * mask_bottom                         &
     213                                ) * ddzw(k) * drho_air(k) * flag
    227214             ENDDO
    228215
    229216!
    230 !--          Vertical diffusion at the first grid point above the surface,
    231 !--          if the momentum flux at the bottom is given by the Prandtl law or
    232 !--          if it is prescribed by the user.
    233 !--          Difference quotient of the momentum flux is not formed over half
    234 !--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
    235 !--          with other (LES) models showed that the values of the momentum
    236 !--          flux becomes too large in this case.
    237 !--          The term containing w(k-1,..) (see above equation) is removed here
    238 !--          because the vertical velocity is assumed to be zero at the surface.
     217!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
     218!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
     219!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
     220!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
     221!--          values of the momentum flux becomes too large in this case.
     222!--          The term containing w(k-1,..) (see above equation) is removed here because the vertical
     223!--          velocity is assumed to be zero at the surface.
    239224             IF ( use_surface_fluxes )  THEN
    240225!
    … …  
    246231                   k   = surf_def_h(0)%k(m)
    247232
    248                    tend(k,j,i) = tend(k,j,i)                                   &
    249                         + ( - ( - surf_def_h(0)%usws(m) )                      &
    250                           ) * ddzw(k) * drho_air(k)
     233                   tend(k,j,i) = tend(k,j,i)                                                       &
     234                                 + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
    251235                ENDDO
    252236!
    … …  
    258242                   k   = surf_def_h(1)%k(m)
    259243
    260                    tend(k,j,i) = tend(k,j,i)                                   &
    261                         + ( - surf_def_h(1)%usws(m)                            &
    262                           ) * ddzw(k) * drho_air(k)
     244                   tend(k,j,i) = tend(k,j,i)                                                       &
     245                                 + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
    263246                ENDDO
    264247!
    … …  
    270253                   k   = surf_lsm_h%k(m)
    271254
    272                    tend(k,j,i) = tend(k,j,i)                                   &
    273                         + ( - ( - surf_lsm_h%usws(m) )                         &
    274                           ) * ddzw(k) * drho_air(k)
     255                   tend(k,j,i) = tend(k,j,i)                                                       &
     256                                 + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
    275257                ENDDO
    276258!
    … …  
    282264                   k   = surf_usm_h%k(m)
    283265
    284                    tend(k,j,i) = tend(k,j,i)                                   &
    285                        + ( - ( - surf_usm_h%usws(m) )                          &
    286                          ) * ddzw(k) * drho_air(k)
     266                   tend(k,j,i) = tend(k,j,i)                                                       &
     267                                 + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
    287268                ENDDO
    288269
    … …  
    297278                   k   = surf_def_h(2)%k(m)
    298279
    299                    tend(k,j,i) = tend(k,j,i)                                   &
    300                         + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
     280                   tend(k,j,i) = tend(k,j,i)                                                       &
     281                                 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
    301282                ENDDO
    302283             ENDIF
    … …  
    308289
    309290
    310 !------------------------------------------------------------------------------!
     291!--------------------------------------------------------------------------------------------------!
    311292! Description:
    312293! ------------
    313294!> Call for grid point i,j
    314 !------------------------------------------------------------------------------!
     295!--------------------------------------------------------------------------------------------------!
    315296    SUBROUTINE diffusion_u_ij( i, j )
    316297
    317        USE arrays_3d,                                                          &
    318            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    319        
    320        USE control_parameters,                                                 &
    321            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    322                   use_top_fluxes
    323        
    324        USE grid_variables,                                                     &
     298       USE arrays_3d,                                                                              &
     299           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
     300
     301       USE control_parameters,                                                                     &
     302           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
     303
     304       USE grid_variables,                                                                         &
    325305           ONLY:  ddx, ddx2, ddy
    326        
    327        USE indices,                                                            &
     306
     307       USE indices,                                                                                &
    328308           ONLY:  nzb, nzt, wall_flags_total_0
    329      
     309
    330310       USE kinds
    331311
    332        USE surface_mod,                                                        &
    333            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    334                    surf_usm_v
     312       USE surface_mod,                                                                            &
     313           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    335314
    336315       IMPLICIT NONE
    … …  
    349328       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
    350329       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
    351        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
    352        REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point 
    353        REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point 
    354        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface 
    355 ! 
     330       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     331       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
     332       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
     333       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
     334!
    356335!--    Compute horizontal diffusion
    357336       DO  k = nzb+1, nzt
    358337!
    359 !--       Predetermine flag to mask topography and wall-bounded grid points. 
    360 !--       It is sufficient to masked only north- and south-facing surfaces, which
    361 !--       need special treatment for the u-component.
    362           flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) ) 
     338!--       Predetermine flag to mask topography and wall-bounded grid points.
     339!--       It is sufficient to masked only north- and south-facing surfaces, which need special
     340!--       treatment for the u-component.
     341          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
    363342          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
    364343          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
    … …  
    368347          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
    369348
    370           tend(k,j,i) = tend(k,j,i)                                            &
    371                        + 2.0_wp * (                                            &
    372                                  km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )     &
    373                                - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )     &
    374                                    ) * ddx2 * flag                             &
    375                                  + (                                           &
    376                   mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy    &
    377                                       + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx    &
    378                                       )                                        &
    379                 - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy    &
    380                                       + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx    &
    381                                       )                                        &
    382                                    ) * ddy  * flag
    383        ENDDO
    384 
    385 !
    386 !--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
    387 !--    surfaces. Note, in the the flat case, loops won't be entered as
    388 !--    start_index > end_index. Furtermore, note, no vertical natural surfaces
    389 !--    so far.           
     349          tend(k,j,i) = tend(k,j,i)                                                                &
     350                                + 2.0_wp * (                                                       &
     351                                         km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                 &
     352                                       - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                 &
     353                                           ) * ddx2 * flag                                         &
     354                                         + (                                                       &
     355                          mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy                &
     356                                              + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                &
     357                                              )                                                    &
     358                        - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy                &
     359                                              + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx                &
     360                                              )                                                    &
     361                                           ) * ddy  * flag
     362       ENDDO
     363
     364!
     365!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in
     366!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
     367!--    vertical natural surfaces so far.
    390368!--    Default-type surfaces
    391369       DO  l = 0, 1
    … …  
    395373             k           = surf_def_v(l)%k(m)
    396374             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
    397           ENDDO   
     375          ENDDO
    398376       ENDDO
    399377!
    … …  
    405383             k           = surf_lsm_v(l)%k(m)
    406384             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
    407           ENDDO   
     385          ENDDO
    408386       ENDDO
    409387!
    … …  
    415393             k           = surf_usm_v(l)%k(m)
    416394             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
    417           ENDDO   
    418        ENDDO
    419 !
    420 !--    Compute vertical diffusion. In case of simulating a surface layer,
    421 !--    respective grid diffusive fluxes are masked (flag 8) within this
    422 !--    loop, and added further below, else, simple gradient approach is
    423 !--    applied. Model top is also mask if top-momentum flux is given.
     395          ENDDO
     396       ENDDO
     397!
     398!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
     399!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
     400!--    approach is applied. Model top is also mask if top-momentum flux is given.
    424401       DO  k = nzb+1, nzt
    425402!
    426 !--       Determine flags to mask topography below and above. Flag 1 is
    427 !--       used to mask topography in general, and flag 8 implies
    428 !--       information about use_surface_fluxes. Flag 9 is used to control
    429 !--       momentum flux at model top.
    430           mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
    431                                BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    432           mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
    433                                BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *     &
    434                         MERGE( 1.0_wp, 0.0_wp,                                 &
    435                                BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    436           flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
    437                                BTEST( wall_flags_total_0(k,j,i), 1 ) )
     403!--       Determine flags to mask topography below and above. Flag 1 is used to mask topography in
     404!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
     405!--       control momentum flux at model top.
     406          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     407          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
     408                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     409          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    438410!
    439411!--       Interpolate eddy diffusivities on staggered gridpoints
    … …  
    441413          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
    442414
    443           tend(k,j,i) = tend(k,j,i)                                            &
    444                         + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
    445                                    + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
    446                                    ) * rho_air_zw(k)   * mask_top              &
    447                           - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
    448                                    + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
    449                                    ) * rho_air_zw(k-1) * mask_bottom           &
     415          tend(k,j,i) = tend(k,j,i)                                                                &
     416                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                       &
     417                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                             &
     418                                   ) * rho_air_zw(k)   * mask_top                                  &
     419                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                       &
     420                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                           &
     421                                   ) * rho_air_zw(k-1) * mask_bottom                               &
    450422                          ) * ddzw(k) * drho_air(k) * flag
    451423       ENDDO
    452424
    453425!
    454 !--    Vertical diffusion at the first surface grid points, if the
    455 !--    momentum flux at the bottom is given by the Prandtl law or if it is
    456 !--    prescribed by the user.
    457 !--    Difference quotient of the momentum flux is not formed over half of
    458 !--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
    459 !--    other (LES) models showed that the values of the momentum flux becomes
    460 !--    too large in this case.
     426!--    Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is
     427!--    given by the Prandtl law or if it is prescribed by the user.
     428!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
     429!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
     430!--    of the momentum flux becomes too large in this case.
    461431       IF ( use_surface_fluxes )  THEN
    462432!
    … …  
    468438             k   = surf_def_h(0)%k(m)
    469439
    470              tend(k,j,i) = tend(k,j,i)                                         &
    471                         + ( - ( - surf_def_h(0)%usws(m) )                      &
    472                           ) * ddzw(k) * drho_air(k)
     440             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
    473441          ENDDO
    474442!
    … …  
    480448             k   = surf_def_h(1)%k(m)
    481449
    482              tend(k,j,i) = tend(k,j,i)                                         &
    483                         + ( - surf_def_h(1)%usws(m)                            &
    484                           ) * ddzw(k) * drho_air(k)
     450             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
    485451          ENDDO
    486452!
    … …  
    492458             k   = surf_lsm_h%k(m)
    493459
    494              tend(k,j,i) = tend(k,j,i)                                         &
    495                         + ( - ( - surf_lsm_h%usws(m) )                         &
    496                           ) * ddzw(k) * drho_air(k)
     460             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
    497461          ENDDO
    498462!
    … …  
    504468             k   = surf_usm_h%k(m)
    505469
    506              tend(k,j,i) = tend(k,j,i)                                         &
    507                        + ( - ( - surf_usm_h%usws(m) )                          &
    508                          ) * ddzw(k) * drho_air(k)
     470             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%usws(m) ) ) * ddzw(k) * drho_air(k)
    509471          ENDDO
    510472
    … …  
    519481             k   = surf_def_h(2)%k(m)
    520482
    521              tend(k,j,i) = tend(k,j,i)                                         &
    522                            + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
     483             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
    523484          ENDDO
    524485       ENDIF

  • palm/trunk/SOURCE/diffusion_v.f90

    r4360 r4583  
    11!> @file diffusion_v.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
    27 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    28 ! topography information used in wall_flags_static_0
    29 !
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
     29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     30! information used in wall_flags_static_0
     31!
    3032! 4329 2019-12-10 15:46:36Z motisi
    3133! Renamed wall_flags_0 to wall_flags_static_0
    32 ! 
     34!
    3335! 4182 2019-08-22 15:20:23Z scharf
    3436! Corrected "Former revisions" section
    35 ! 
     37!
    3638! 3655 2019-01-07 16:51:22Z knoop
    3739! OpenACC port for SPEC
    … …  
    4446! ------------
    4547!> Diffusion term of the v-component
    46 !------------------------------------------------------------------------------!
     48!--------------------------------------------------------------------------------------------------!
    4749 MODULE diffusion_v_mod
    48  
     50
    4951
    5052    PRIVATE
    … …  
    5961
    6062
    61 !------------------------------------------------------------------------------!
     63!--------------------------------------------------------------------------------------------------!
    6264! Description:
    6365! ------------
    6466!> Call for all grid points
    65 !------------------------------------------------------------------------------!
     67!--------------------------------------------------------------------------------------------------!
    6668    SUBROUTINE diffusion_v
    6769
    68        USE arrays_3d,                                                          &
    69            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    70        
    71        USE control_parameters,                                                 &
    72            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    73                   use_top_fluxes
    74        
    75        USE grid_variables,                                                     &
     70       USE arrays_3d,                                                                              &
     71           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
     72
     73       USE control_parameters,                                                                     &
     74           ONLY:  constant_top_momentumflux, use_surface_fluxes,  use_top_fluxes
     75
     76       USE grid_variables,                                                                         &
    7677           ONLY:  ddx, ddy, ddy2
    77        
    78        USE indices,                                                            &
     78
     79       USE indices,                                                                                &
    7980           ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
    80        
     81
    8182       USE kinds
    8283
    83        USE surface_mod,                                                        &
    84            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    85                    surf_usm_v
     84       USE surface_mod,                                                                            &
     85           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    8686
    8787       IMPLICIT NONE
    … …  
    100100       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
    101101       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
    102        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
    103        REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
    104        REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
    105        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface     
     102       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     103       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     104       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
     105       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    106106
    107107       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
    … …  
    122122
    123123!
    124 !--             Predetermine flag to mask topography and wall-bounded grid points. 
    125 !--             It is sufficient to masked only east- and west-facing surfaces, which
    126 !--             need special treatment for the v-component.
    127                 flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
     124!--             Predetermine flag to mask topography and wall-bounded grid points.
     125!--             It is sufficient to masked only east- and west-facing surfaces, which need special
     126!--             treatment for the v-component.
     127                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
    128128                mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
    129129                mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
    … …  
    133133                kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    134134
    135                 tend(k,j,i) = tend(k,j,i) +    (                             &
    136                           mask_east * kmxp * (                               &
    137                                  ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    138                                + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    139                                              )                               &
    140                         - mask_west * kmxm * (                               &
    141                                  ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    142                                + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    143                                              )                               &
    144                                                ) * ddx  * flag               &
    145                                     + 2.0_wp * (                             &
    146                                   km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
    147                                 - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
    148                                                ) * ddy2 * flag
     135                tend(k,j,i) = tend(k,j,i) +          (                                             &
     136                                mask_east * kmxp * (                                               &
     137                                       ( v(k,j,i+1) - v(k,j,i)     ) * ddx                         &
     138                                     + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                         &
     139                                                   )                                               &
     140                              - mask_west * kmxm * (                                               &
     141                                       ( v(k,j,i) - v(k,j,i-1) ) * ddx                             &
     142                                     + ( u(k,j,i) - u(k,j-1,i) ) * ddy                             &
     143                                                   )                                               &
     144                                                     ) * ddx  * flag                               &
     145                                          + 2.0_wp * (                                             &
     146                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                  &
     147                                      - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                  &
     148                                                     ) * ddy2 * flag
    149149
    150150             ENDDO
    151151
    152152!
    153 !--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
    154 !--          surfaces. Note, in the the flat case, loops won't be entered as
    155 !--          start_index > end_index. Furtermore, note, no vertical natural surfaces
    156 !--          so far.           
     153!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note,
     154!--          in the the flat case, loops won't be entered as start_index > end_index. Furtermore,
     155!--          note, no vertical natural surfaces so far.
    157156!--          Default-type surfaces
    158157             DO  l = 2, 3
    … …  
    161160                DO  m = surf_s, surf_e
    162161                   k           = surf_def_v(l)%k(m)
    163                    tend(k,j,i) = tend(k,j,i) +                                 &
    164                                     surf_def_v(l)%mom_flux_uv(m) * ddx
    165                 ENDDO   
     162                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
     163                ENDDO
    166164             ENDDO
    167165!
    … …  
    172170                DO  m = surf_s, surf_e
    173171                   k           = surf_lsm_v(l)%k(m)
    174                    tend(k,j,i) = tend(k,j,i) +                                 &
    175                                     surf_lsm_v(l)%mom_flux_uv(m) * ddx
    176                 ENDDO   
     172                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
     173                ENDDO
    177174             ENDDO
    178175!
    … …  
    183180                DO  m = surf_s, surf_e
    184181                   k           = surf_usm_v(l)%k(m)
    185                    tend(k,j,i) = tend(k,j,i) +                                 &
    186                                     surf_usm_v(l)%mom_flux_uv(m) * ddx
    187                 ENDDO   
     182                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
     183                ENDDO
    188184             ENDDO
    189185!
    190 !--          Compute vertical diffusion. In case of simulating a surface layer,
    191 !--          respective grid diffusive fluxes are masked (flag 10) within this
    192 !--          loop, and added further below, else, simple gradient approach is
    193 !--          applied. Model top is also mask if top-momentum flux is given.
     186!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
     187!--          diffusive fluxes are masked (flag 10) within this loop, and added further below, else,
     188!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
     189!--          given.
    194190             DO  k = nzb+1, nzt
    195191!
    196 !--             Determine flags to mask topography below and above. Flag 2 is
    197 !--             used to mask topography in general, while flag 8 implies also
    198 !--             information about use_surface_fluxes. Flag 9 is used to control
    199 !--             momentum flux at model top. 
    200                 mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
    201                                  BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    202                 mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
    203                                  BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
    204                               MERGE( 1.0_wp, 0.0_wp,                           &
    205                                  BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    206                 flag        = MERGE( 1.0_wp, 0.0_wp,                           &
    207                                  BTEST( wall_flags_total_0(k,j,i), 2 ) )
     192!--             Determine flags to mask topography below and above. Flag 2 is used to mask
     193!--             topography in general, while flag 8 implies also information about
     194!--             use_surface_fluxes. Flag 9 is used to control momentum flux at model top.
     195                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     196                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
     197                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     198                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    208199!
    209200!--             Interpolate eddy diffusivities on staggered gridpoints
    210                 kmzp = 0.25_wp * &
    211                        ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    212                 kmzm = 0.25_wp * &
    213                        ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    214 
    215                 tend(k,j,i) = tend(k,j,i)                                      &
    216                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    217                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    218                       &            ) * rho_air_zw(k)   * mask_top              &
    219                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    220                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    221                       &            ) * rho_air_zw(k-1) * mask_bottom           &
    222                       &   ) * ddzw(k) * drho_air(k) * flag
     201                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
     202                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
     203
     204                tend(k,j,i) = tend(k,j,i)                                                          &
     205                              & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                 &
     206                              &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                       &
     207                              &            ) * rho_air_zw(k)   * mask_top                          &
     208                              &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)               &
     209                              &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                   &
     210                              &            ) * rho_air_zw(k-1) * mask_bottom                       &
     211                              &   ) * ddzw(k) * drho_air(k) * flag
    223212             ENDDO
    224213
    225214!
    226 !--          Vertical diffusion at the first grid point above the surface,
    227 !--          if the momentum flux at the bottom is given by the Prandtl law
    228 !--          or if it is prescribed by the user.
    229 !--          Difference quotient of the momentum flux is not formed over
    230 !--          half of the grid spacing (2.0*ddzw(k)) any more, since the
    231 !--          comparison with other (LES) models showed that the values of
    232 !--          the momentum flux becomes too large in this case.
     215!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
     216!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
     217!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
     218!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
     219!--          values of the momentum flux becomes too large in this case.
    233220             IF ( use_surface_fluxes )  THEN
    234221!
    … …  
    239226                   k   = surf_def_h(0)%k(m)
    240227
    241                    tend(k,j,i) = tend(k,j,i)                                   &
    242                         + ( - ( - surf_def_h(0)%vsws(m) )                      &
    243                           ) * ddzw(k) * drho_air(k)
     228                   tend(k,j,i) = tend(k,j,i)                                                       &
     229                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
    244230                ENDDO
    245231!
    … …  
    250236                   k   = surf_def_h(1)%k(m)
    251237
    252                    tend(k,j,i) = tend(k,j,i)                                   &
    253                         + ( - surf_def_h(1)%vsws(m)                            &
    254                           ) * ddzw(k) * drho_air(k)
     238                   tend(k,j,i) = tend(k,j,i)                                                       &
     239                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
    255240                ENDDO
    256241!
    … …  
    261246                   k   = surf_lsm_h%k(m)
    262247
    263                    tend(k,j,i) = tend(k,j,i)                                   &
    264                         + ( - ( - surf_lsm_h%vsws(m) )                         &
    265                           ) * ddzw(k) * drho_air(k)
     248                   tend(k,j,i) = tend(k,j,i)                                                       &
     249                                 + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    266250
    267251                ENDDO
    … …  
    273257                   k   = surf_usm_h%k(m)
    274258
    275                    tend(k,j,i) = tend(k,j,i)                                   &
    276                         + ( - ( - surf_usm_h%vsws(m) )                         &
    277                           ) * ddzw(k) * drho_air(k)
     259                   tend(k,j,i) = tend(k,j,i)                                                       &
     260                                 + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    278261
    279262                ENDDO
    … …  
    288271                   k   = surf_def_h(2)%k(m)
    289272
    290                    tend(k,j,i) = tend(k,j,i)                                   &
    291                            + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     273                   tend(k,j,i) = tend(k,j,i)                                                       &
     274                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
    292275                ENDDO
    293276             ENDIF
    … …  
    299282
    300283
    301 !------------------------------------------------------------------------------!
     284!--------------------------------------------------------------------------------------------------!
    302285! Description:
    303286! ------------
    304287!> Call for grid point i,j
    305 !------------------------------------------------------------------------------!
     288!--------------------------------------------------------------------------------------------------!
    306289    SUBROUTINE diffusion_v_ij( i, j )
    307290
    308        USE arrays_3d,                                                          &
    309            ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
    310        
    311        USE control_parameters,                                                 &
    312            ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
    313                   use_top_fluxes
    314        
    315        USE grid_variables,                                                     &
     291       USE arrays_3d,                                                                              &
     292           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
     293
     294       USE control_parameters,                                                                     &
     295           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
     296
     297       USE grid_variables,                                                                         &
    316298           ONLY:  ddx, ddy, ddy2
    317        
    318        USE indices,                                                            &
     299
     300       USE indices,                                                                                &
    319301           ONLY:  nzb, nzt, wall_flags_total_0
    320        
     302
    321303       USE kinds
    322304
    323        USE surface_mod,                                                        &
    324            ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
    325                    surf_usm_v
     305       USE surface_mod,                                                                            &
     306           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
    326307
    327308       IMPLICIT NONE
    … …  
    341322       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
    342323       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
    343        REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface 
    344        REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point 
    345        REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point 
     324       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
     325       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
     326       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
    346327       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
    347328
    … …  
    350331       DO  k = nzb+1, nzt
    351332!
    352 !--       Predetermine flag to mask topography and wall-bounded grid points. 
    353 !--       It is sufficient to masked only east- and west-facing surfaces, which
    354 !--       need special treatment for the v-component.
    355           flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) ) 
     333!--       Predetermine flag to mask topography and wall-bounded grid points.
     334!--       It is sufficient to masked only east- and west-facing surfaces, which need special
     335!--       treatment for the v-component.
     336          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
    356337          mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
    357338          mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
    … …  
    361342          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
    362343
    363           tend(k,j,i) = tend(k,j,i) +          (                             &
    364                           mask_east * kmxp * (                               &
    365                                  ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
    366                                + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
    367                                              )                               &
    368                         - mask_west * kmxm * (                               &
    369                                  ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
    370                                + ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
    371                                              )                               &
    372                                                ) * ddx  * flag               &
    373                                     + 2.0_wp * (                             &
    374                                   km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )  &
    375                                 - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )  &
     344          tend(k,j,i) = tend(k,j,i) +          (                                                   &
     345                          mask_east * kmxp * (                                                     &
     346                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx                               &
     347                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                               &
     348                                             )                                                     &
     349                        - mask_west * kmxm * (                                                     &
     350                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx                                   &
     351                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy                                   &
     352                                             )                                                     &
     353                                               ) * ddx  * flag                                     &
     354                                    + 2.0_wp * (                                                   &
     355                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                        &
     356                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                        &
    376357                                               ) * ddy2 * flag
    377358       ENDDO
    378359
    379360!
    380 !--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3)
    381 !--    surfaces. Note, in the the flat case, loops won't be entered as
    382 !--    start_index > end_index. Furtermore, note, no vertical natural surfaces
    383 !--    so far.           
     361!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note, in the
     362!--    the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
     363!--    vertical natural surfaces so far.
    384364!--    Default-type surfaces
    385365       DO  l = 2, 3
    … …  
    389369             k           = surf_def_v(l)%k(m)
    390370             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
    391           ENDDO   
     371          ENDDO
    392372       ENDDO
    393373!
    … …  
    399379             k           = surf_lsm_v(l)%k(m)
    400380             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
    401           ENDDO   
     381          ENDDO
    402382       ENDDO
    403383!
    … …  
    409389             k           = surf_usm_v(l)%k(m)
    410390             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
    411           ENDDO   
    412        ENDDO
    413 !
    414 !--    Compute vertical diffusion. In case of simulating a surface layer,
    415 !--    respective grid diffusive fluxes are masked (flag 8) within this
    416 !--    loop, and added further below, else, simple gradient approach is
    417 !--    applied. Model top is also mask if top-momentum flux is given.
     391          ENDDO
     392       ENDDO
     393!
     394!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
     395!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
     396!--    approach is applied. Model top is also mask if top-momentum flux is given.
    418397       DO  k = nzb+1, nzt
    419398!
    420 !--       Determine flags to mask topography below and above. Flag 2 is
    421 !--       used to mask topography in general, while flag 10 implies also
    422 !--       information about use_surface_fluxes. Flag 9 is used to control
    423 !--       momentum flux at model top. 
    424           mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
    425                                BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
    426           mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
    427                                BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *     &
    428                         MERGE( 1.0_wp, 0.0_wp,                                 &
    429                                BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
    430           flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
    431                                BTEST( wall_flags_total_0(k,j,i), 2 ) )
     399!--       Determine flags to mask topography below and above. Flag 2 is used to mask topography in
     400!--       general, while flag 10 implies also information about use_surface_fluxes. Flag 9 is used
     401!--       to control momentum flux at model top.
     402          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
     403          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
     404                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
     405          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    432406!
    433407!--       Interpolate eddy diffusivities on staggered gridpoints
    … …  
    435409          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
    436410
    437           tend(k,j,i) = tend(k,j,i)                                            &
    438                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
    439                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
    440                       &            ) * rho_air_zw(k)   * mask_top              &
    441                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
    442                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
    443                       &            ) * rho_air_zw(k-1) * mask_bottom           &
    444                       &   ) * ddzw(k) * drho_air(k) * flag
    445        ENDDO
    446 
    447 !
    448 !--    Vertical diffusion at the first grid point above the surface, if the
    449 !--    momentum flux at the bottom is given by the Prandtl law or if it is
    450 !--    prescribed by the user.
    451 !--    Difference quotient of the momentum flux is not formed over half of
    452 !--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
    453 !--    other (LES) models showed that the values of the momentum flux becomes
    454 !--    too large in this case.
     411          tend(k,j,i) = tend(k,j,i)                                                                &
     412                        & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                       &
     413                        &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                             &
     414                        &            ) * rho_air_zw(k)   * mask_top                                &
     415                        &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)                     &
     416                        &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                         &
     417                        &            ) * rho_air_zw(k-1) * mask_bottom                             &
     418                        &   ) * ddzw(k) * drho_air(k) * flag
     419       ENDDO
     420
     421!
     422!--    Vertical diffusion at the first grid point above the surface, if the momentum flux at the
     423!--    bottom is given by the Prandtl law or if it is prescribed by the user.
     424!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
     425!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
     426!--    of the momentum flux becomes too large in this case.
    455427       IF ( use_surface_fluxes )  THEN
    456428!
    … …  
    461433             k   = surf_def_h(0)%k(m)
    462434
    463              tend(k,j,i) = tend(k,j,i)                                         &
    464                         + ( - ( - surf_def_h(0)%vsws(m) )                      &
    465                           ) * ddzw(k) * drho_air(k)
     435             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
    466436          ENDDO
    467437!
    … …  
    472442             k   = surf_def_h(1)%k(m)
    473443
    474              tend(k,j,i) = tend(k,j,i)                                         &
    475                         + ( - surf_def_h(1)%vsws(m)                            &
    476                           ) * ddzw(k) * drho_air(k)
     444             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
    477445          ENDDO
    478446!
    … …  
    483451             k   = surf_lsm_h%k(m)
    484452
    485              tend(k,j,i) = tend(k,j,i)                                         &
    486                         + ( - ( - surf_lsm_h%vsws(m) )                         &
    487                           ) * ddzw(k) * drho_air(k)
     453             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    488454
    489455          ENDDO
    … …  
    495461             k   = surf_usm_h%k(m)
    496462
    497              tend(k,j,i) = tend(k,j,i)                                         &
    498                         + ( - ( - surf_usm_h%vsws(m) )                         &
    499                           ) * ddzw(k) * drho_air(k)
     463             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
    500464
    501465          ENDDO
    … …  
    510474             k   = surf_def_h(2)%k(m)
    511475
    512              tend(k,j,i) = tend(k,j,i)                                         &
    513                            + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
     476             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
    514477          ENDDO
    515478       ENDIF

  • palm/trunk/SOURCE/diffusion_w.f90

    r4360 r4583  
    11!> @file diffusion_w.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
    27 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    28 ! topography information used in wall_flags_static_0
    29 !
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
     29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     30! information used in wall_flags_static_0
     31!
    3032! 4329 2019-12-10 15:46:36Z motisi
    3133! Renamed wall_flags_0 to wall_flags_static_0
    32 ! 
     34!
    3335! 4182 2019-08-22 15:20:23Z scharf
    3436! Corrected "Former revisions" section
    35 ! 
     37!
    3638! 3655 2019-01-07 16:51:22Z knoop
    3739! OpenACC port for SPEC
    … …  
    4446! ------------
    4547!> Diffusion term of the w-component
    46 !------------------------------------------------------------------------------!
     48!--------------------------------------------------------------------------------------------------!
    4749 MODULE diffusion_w_mod
    48  
     50
    4951
    5052    PRIVATE
    … …  
    5961
    6062
    61 !------------------------------------------------------------------------------!
     63!--------------------------------------------------------------------------------------------------!
    6264! Description:
    6365! ------------
    6466!> Call for all grid points
    65 !------------------------------------------------------------------------------!
     67!--------------------------------------------------------------------------------------------------!
    6668    SUBROUTINE diffusion_w
    6769
    68        USE arrays_3d,                                                          &         
    69            ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    70            
    71        USE grid_variables,                                                     &     
     70       USE arrays_3d,                                                                              &
     71           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
     72
     73       USE grid_variables,                                                                         &
    7274           ONLY :  ddx, ddy
    73            
    74        USE indices,                                                            &           
     75
     76       USE indices,                                                                                &
    7577           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
    76            
     78
    7779       USE kinds
    7880
    79        USE surface_mod,                                                        &
     81       USE surface_mod,                                                                            &
    8082           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
    8183
    … …  
    8991       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
    9092       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
    91        
     93
    9294       REAL(wp) ::  flag              !< flag to mask topography grid points
    9395       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
    … …  
    9597       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
    9698       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
     99       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
     100       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
     101       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
    97102       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
    98        REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
    99        REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
    100        REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
    101103
    102104
    … …  
    116118             DO  k = nzb+1, nzt-1
    117119!
    118 !--             Predetermine flag to mask topography and wall-bounded grid points.
    119                 flag       = MERGE( 1.0_wp, 0.0_wp,                            &
    120                                     BTEST( wall_flags_total_0(k,j,i),   3 ) )
    121                 mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
    122                                     BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
    123                 mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
    124                                     BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
    125                 mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
    126                                     BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
    127                 mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
    128                                     BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
     120!--             Predetermine flag to mask topography and wall-bounded grid points.
     121                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
     122                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
     123                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
     124                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
     125                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
    129126!
    130127!--             Interpolate eddy diffusivities on staggered gridpoints
    131                 kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +               &
     128                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +                                   &
    132129                                   km(k+1,j,i) + km(k+1,j,i+1) )
    133                 kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +               &
     130                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +                                   &
    134131                                   km(k+1,j,i) + km(k+1,j,i-1) )
    135                 kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
     132                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
    136133                                   km(k,j+1,i) + km(k+1,j+1,i) )
    137                 kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
     134                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
    138135                                   km(k,j-1,i) + km(k+1,j-1,i) )
    139136
    140                 tend(k,j,i) = tend(k,j,i)                                      &
    141                        + ( mask_east *  kmxp * (                               &
    142                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    143                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    144                                                )                               &
    145                          - mask_west * kmxm *  (                               &
    146                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
    147                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
    148                                                )                               &
    149                          ) * ddx                                 * flag        &
    150                        + ( mask_north * kmyp * (                               &
    151                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    152                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    153                                                )                               &
    154                          - mask_south * kmym * (                               &
    155                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
    156                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
    157                                                )                               &
    158                          ) * ddy                                 * flag        &
    159                        + 2.0_wp * (                                            &
    160                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
    161                                      * rho_air(k+1)                            &
    162                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
    163                                      * rho_air(k)                              &
    164                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
    165              ENDDO
    166 
    167 !
    168 !--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
    169 !--          surfaces. Note, in the the flat case, loops won't be entered as
    170 !--          start_index > end_index. Furtermore, note, no vertical natural surfaces
    171 !--          so far.           
     137                tend(k,j,i) = tend(k,j,i)                                                          &
     138                              + ( mask_east *  kmxp * (                                            &
     139                                          ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                      &
     140                                        + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                &
     141                                                      )                                            &
     142                                - mask_west * kmxm *  (                                            &
     143                                          ( w(k,j,i)     - w(k,j,i-1) ) * ddx                      &
     144                                        + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                &
     145                                                      )                                            &
     146                                ) * ddx                                 * flag                     &
     147                              + ( mask_north * kmyp * (                                            &
     148                                          ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                      &
     149                                        + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                &
     150                                                      )                                            &
     151                                - mask_south * kmym * (                                            &
     152                                          ( w(k,j,i)     - w(k,j-1,i) ) * ddy                      &
     153                                        + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                &
     154                                                      )                                            &
     155                                ) * ddy                                 * flag                     &
     156                              + 2.0_wp * (                                                         &
     157                                km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1)              &
     158                                            * rho_air(k+1)                                         &
     159                              - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                &
     160                                            * rho_air(k)                                           &
     161                                         ) * ddzu(k+1) * drho_air_zw(k) * flag
     162             ENDDO
     163
     164!
     165!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces.
     166!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
     167!--          Furtermore, note, no vertical natural surfaces so far.
    172168!--          Default-type surfaces
    173169             DO  l = 0, 1
    … …  
    176172                DO  m = surf_s, surf_e
    177173                   k           = surf_def_v(l)%k(m)
    178                    tend(k,j,i) = tend(k,j,i) +                                 &
    179                                      surf_def_v(l)%mom_flux_w(m) * ddy
    180                 ENDDO   
     174                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
     175                ENDDO
    181176             ENDDO
    182177!
    … …  
    187182                DO  m = surf_s, surf_e
    188183                   k           = surf_lsm_v(l)%k(m)
    189                    tend(k,j,i) = tend(k,j,i) +                                 &
    190                                      surf_lsm_v(l)%mom_flux_w(m) * ddy
    191                 ENDDO   
     184                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
     185                ENDDO
    192186             ENDDO
    193187!
    … …  
    198192                DO  m = surf_s, surf_e
    199193                   k           = surf_usm_v(l)%k(m)
    200                    tend(k,j,i) = tend(k,j,i) +                                 &
    201                                      surf_usm_v(l)%mom_flux_w(m) * ddy
    202                 ENDDO   
    203              ENDDO
    204 !
    205 !--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
    206 !--          surface.
     194                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
     195                ENDDO
     196             ENDDO
     197!
     198!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface.
    207199!--          Default-type surfaces
    208200             DO  l = 2, 3
    … …  
    211203                DO  m = surf_s, surf_e
    212204                   k           = surf_def_v(l)%k(m)
    213                    tend(k,j,i) = tend(k,j,i) +                                 &
    214                                      surf_def_v(l)%mom_flux_w(m) * ddx
    215                 ENDDO   
     205                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
     206                ENDDO
    216207             ENDDO
    217208!
    … …  
    222213                DO  m = surf_s, surf_e
    223214                   k           = surf_lsm_v(l)%k(m)
    224                    tend(k,j,i) = tend(k,j,i) +                                 &
    225                                      surf_lsm_v(l)%mom_flux_w(m) * ddx
    226                 ENDDO   
     215                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
     216                ENDDO
    227217             ENDDO
    228218!
    … …  
    233223                DO  m = surf_s, surf_e
    234224                   k           = surf_usm_v(l)%k(m)
    235                    tend(k,j,i) = tend(k,j,i) +                                 &
    236                                      surf_usm_v(l)%mom_flux_w(m) * ddx
    237                 ENDDO   
     225                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
     226                ENDDO
    238227             ENDDO
    239228
    … …  
    244233
    245234
    246 !------------------------------------------------------------------------------!
     235!--------------------------------------------------------------------------------------------------!
    247236! Description:
    248237! ------------
    249238!> Call for grid point i,j
    250 !------------------------------------------------------------------------------!
     239!--------------------------------------------------------------------------------------------------!
    251240    SUBROUTINE diffusion_w_ij( i, j )
    252241
    253        USE arrays_3d,                                                          &         
    254            ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
    255            
    256        USE grid_variables,                                                     &     
     242       USE arrays_3d,                                                                              &
     243           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
     244
     245       USE grid_variables,                                                                         &
    257246           ONLY :  ddx, ddy
    258            
    259        USE indices,                                                            &           
     247
     248       USE indices,                                                                                &
    260249           ONLY :  nzb, nzt, wall_flags_total_0
    261            
     250
    262251       USE kinds
    263252
    264        USE surface_mod,                                                        &
     253       USE surface_mod,                                                                            &
    265254           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
    266255
    … …  
    275264       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
    276265       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
    277        
     266
    278267       REAL(wp) ::  flag              !< flag to mask topography grid points
    279268       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
    … …  
    281270       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
    282271       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
     272       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
     273       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
     274       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
    283275       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
    284        REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
    285        REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
    286        REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
    287276
    288277
    289278       DO  k = nzb+1, nzt-1
    290279!
    291 !--       Predetermine flag to mask topography and wall-bounded grid points. 
    292           flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) ) 
     280!--       Predetermine flag to mask topography and wall-bounded grid points.
     281          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
    293282          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
    294283          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
    … …  
    302291          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
    303292
    304           tend(k,j,i) = tend(k,j,i)                                            &
    305                        + ( mask_east *  kmxp * (                               &
    306                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
    307                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
    308                                                )                               &
    309                          - mask_west * kmxm *  (                               &
    310                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
    311                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
    312                                                )                               &
    313                          ) * ddx                                 * flag        &
    314                        + ( mask_north * kmyp * (                               &
    315                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
    316                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
    317                                                )                               &
    318                          - mask_south * kmym * (                               &
    319                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
    320                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
    321                                                )                               &
    322                          ) * ddy                                 * flag        &
    323                        + 2.0_wp * (                                            &
    324                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
    325                                      * rho_air(k+1)                            &
    326                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
    327                                      * rho_air(k)                              &
    328                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
    329        ENDDO
    330 !
    331 !--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
    332 !--    surfaces. Note, in the the flat case, loops won't be entered as
    333 !--    start_index > end_index. Furtermore, note, no vertical natural surfaces
    334 !--    so far.           
     293          tend(k,j,i) = tend(k,j,i)                                                                &
     294                        + ( mask_east *  kmxp * (                                                  &
     295                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                            &
     296                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                      &
     297                                                )                                                  &
     298                          - mask_west * kmxm *  (                                                  &
     299                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx                            &
     300                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                      &
     301                                                )                                                  &
     302                          ) * ddx                                 * flag                           &
     303                        + ( mask_north * kmyp * (                                                  &
     304                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                            &
     305                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                      &
     306                                                )                                                  &
     307                          - mask_south * kmym * (                                                  &
     308                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy                            &
     309                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                      &
     310                                                )                                                  &
     311                          ) * ddy                                 * flag                           &
     312                        + 2.0_wp * (                                                               &
     313                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)                      &
     314                                      * rho_air(k+1)                                               &
     315                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                      &
     316                                      * rho_air(k)                                                 &
     317                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
     318       ENDDO
     319!
     320!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in
     321!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
     322!--    vertical natural surfaces so far.
    335323!--    Default-type surfaces
    336324       DO  l = 0, 1
    … …  
    339327          DO  m = surf_s, surf_e
    340328             k           = surf_def_v(l)%k(m)
    341              tend(k,j,i) = tend(k,j,i) +                                       &
    342                                      surf_def_v(l)%mom_flux_w(m) * ddy
    343           ENDDO   
     329             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
     330          ENDDO
    344331       ENDDO
    345332!
    … …  
    350337          DO  m = surf_s, surf_e
    351338             k           = surf_lsm_v(l)%k(m)
    352              tend(k,j,i) = tend(k,j,i) +                                       &
    353                                      surf_lsm_v(l)%mom_flux_w(m) * ddy
    354           ENDDO   
     339             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
     340          ENDDO
    355341       ENDDO
    356342!
    … …  
    361347          DO  m = surf_s, surf_e
    362348             k           = surf_usm_v(l)%k(m)
    363              tend(k,j,i) = tend(k,j,i) +                                       &
    364                                      surf_usm_v(l)%mom_flux_w(m) * ddy
    365           ENDDO   
    366        ENDDO
    367 !
    368 !--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
    369 !--    surfaces.
     349             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
     350          ENDDO
     351       ENDDO
     352!
     353!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces.
    370354!--    Default-type surfaces
    371355       DO  l = 2, 3
    … …  
    374358          DO  m = surf_s, surf_e
    375359             k           = surf_def_v(l)%k(m)
    376              tend(k,j,i) = tend(k,j,i) +                                       &
    377                                      surf_def_v(l)%mom_flux_w(m) * ddx
    378           ENDDO   
     360             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
     361          ENDDO
    379362       ENDDO
    380363!
    … …  
    385368          DO  m = surf_s, surf_e
    386369             k           = surf_lsm_v(l)%k(m)
    387              tend(k,j,i) = tend(k,j,i) +                                       &
    388                                      surf_lsm_v(l)%mom_flux_w(m) * ddx
    389           ENDDO   
     370             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
     371          ENDDO
    390372       ENDDO
    391373!
    … …  
    396378          DO  m = surf_s, surf_e
    397379             k           = surf_usm_v(l)%k(m)
    398              tend(k,j,i) = tend(k,j,i) +                                       &
    399                                      surf_usm_v(l)%mom_flux_w(m) * ddx
    400           ENDDO   
     380             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
     381          ENDDO
    401382       ENDDO
    402383

  • palm/trunk/SOURCE/disturb_field.f90

    r4457 r4583  
    11!> @file disturb_field.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4457 2020-03-11 14:20:43Z raasch
    2729! use statement for exchange horiz added
    28 ! 
     30!
    2931! 4360 2020-01-07 11:25:50Z suehring
    30 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    31 ! topography information used in wall_flags_static_0
    32 ! 
     32! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     33! information used in wall_flags_static_0
     34!
    3335! 4329 2019-12-10 15:46:36Z motisi
    3436! Renamed wall_flags_0 to wall_flags_static_0
    35 ! 
     37!
    3638! 4237 2019-09-25 11:33:42Z knoop
    3739! Added missing OpenMP directives
    38 ! 
     40!
    3941! 4182 2019-08-22 15:20:23Z scharf
    4042! Corrected "Former revisions" section
    41 ! 
     43!
    4244! 3849 2019-04-01 16:35:16Z knoop
    4345! Corrected "Former revisions" section
    … …  
    5052! ------------
    5153!> Imposing a random perturbation on a 3D-array.
    52 !> On parallel computers, the random number generator is as well called for all
    53 !> gridpoints of the total domain to ensure, regardless of the number of PEs
    54 !> used, that the elements of the array have the same values in the same
    55 !> order in every case. The perturbation range is steered by dist_range.
    56 !------------------------------------------------------------------------------!
     54!> On parallel computers, the random number generator is as well called for all gridpoints of the
     55!> total domain to ensure, regardless of the number of PEs used, that the elements of the array have
     56!> the same values in the same order in every case. The perturbation range is steered by dist_range.
     57!--------------------------------------------------------------------------------------------------!
    5758 SUBROUTINE disturb_field( var_char, dist1, field )
    58  
    59 
    60     USE control_parameters,                                                    &
    61         ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
    62                disturbance_amplitude, disturbance_created,                     &
    63                disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
     59
     60
     61    USE control_parameters,                                                                        &
     62        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range, disturbance_amplitude,          &
     63               disturbance_created, disturbance_level_ind_b, disturbance_level_ind_t, iran,        &
    6464               random_generator, topography
    65                
    66     USE cpulog,                                                                &
     65
     66    USE cpulog,                                                                                    &
    6767        ONLY:  cpu_log, log_point
    68        
    69     USE exchange_horiz_mod,                                                    &
     68
     69    USE exchange_horiz_mod,                                                                        &
    7070        ONLY:  exchange_horiz
    7171
    72     USE indices,                                                               &
    73         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
    74                nzt, wall_flags_total_0
    75        
     72    USE indices,                                                                                   &
     73        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, nzt,                &
     74               wall_flags_total_0
     75
    7676    USE kinds
    77    
    78     USE random_function_mod,                                                   &
     77
     78    USE random_function_mod,                                                                       &
    7979        ONLY: random_function
    80        
    81     USE random_generator_parallel,                                             &
    82         ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
    83                seq_random_array
     80
     81    USE random_generator_parallel,                                                                 &
     82        ONLY:  random_number_parallel, random_seed_parallel, random_dummy, seq_random_array
    8483
    8584    IMPLICIT NONE
    … …  
    8786    CHARACTER (LEN = *) ::  var_char !< flag to distinguish betwenn u- and v-component
    8887
    89     INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid 
     88    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid
    9089    INTEGER(iwp) ::  i       !< index variable
    9190    INTEGER(iwp) ::  j       !< index variable
    … …  
    9392
    9493    REAL(wp) ::  randomnumber  !<
    95    
     94
    9695    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
    9796    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
    98    
     97
    9998    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
    10099
    … …  
    105104    flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' )
    106105!
    107 !-- Create an additional temporary array and initialize the arrays needed
    108 !-- to store the disturbance
     106!-- Create an additional temporary array and initialize the arrays needed to store the disturbance
    109107    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    110108!$ACC DATA CREATE(dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
    … …  
    123121          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
    124122             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
    125                 randomnumber = 3.0_wp * disturbance_amplitude *                &
    126                                ( random_function( iran ) - 0.5_wp )
    127                 IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
    128                      nyn >= j )                                                &
    129                 THEN
     123                randomnumber = 3.0_wp * disturbance_amplitude * ( random_function( iran ) - 0.5_wp )
     124                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND. nyn >= j )  THEN
    130125                   dist1(k,j,i) = randomnumber
    131126                ENDIF
    … …  
    139134             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
    140135                CALL random_number_parallel( random_dummy )
    141                 randomnumber = 3.0_wp * disturbance_amplitude *                &
    142                                ( random_dummy - 0.5_wp )
    143                 IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
    144                      nyn >= j )                                                &
    145                 THEN
     136                randomnumber = 3.0_wp * disturbance_amplitude * ( random_dummy - 0.5_wp )
     137                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
    146138                   dist1(k,j,i) = randomnumber
    147139                ENDIF
    … …  
    155147             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
    156148                CALL RANDOM_NUMBER( randomnumber )
    157                 randomnumber = 3.0_wp * disturbance_amplitude *                &
    158                                 ( randomnumber - 0.5_wp )
    159                 IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
    160                 THEN
     149                randomnumber = 3.0_wp * disturbance_amplitude * ( randomnumber - 0.5_wp )
     150                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
    161151                   dist1(k,j,i) = randomnumber
    162152                ENDIF
    … …  
    177167!
    178168!-- Applying the Shuman filter in order to smooth the perturbations.
    179 !-- Neighboured grid points in all three directions are used for the
    180 !-- filter operation.
    181 !-- Loop has been splitted to make runs reproducible on HLRN systems using
    182 !-- compiler option -O3
     169!-- Neighboured grid points in all three directions are used for the filter operation.
     170!-- Loop has been splitted to make runs reproducible on HLRN systems using compiler option -O3
    183171     !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
    184172     !$OMP PARALLEL DO PRIVATE(i, j, k)
    … …  
    186174        DO  j = nys, nyn
    187175          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
    188              dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
    189                             + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
     176             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                                      &
     177                            + dist1(k,j+1,i) + dist1(k+1,j,i)                                      &
    190178                            ) / 12.0_wp
    191179          ENDDO
    192180          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
    193               dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
    194                             + 6.0_wp * dist1(k,j,i)                            &
    195                             ) / 12.0_wp
     181              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)                      &
     182                            + 6.0_wp * dist1(k,j,i)                                                &
     183                                            ) / 12.0_wp
    196184          ENDDO
    197185        ENDDO
    … …  
    208196       DO  j = nys, nyn
    209197          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
    210              dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
    211                             + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
    212                             + 6.0_wp * dist2(k,j,i)                            &
     198             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i)                     &
     199                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i)                     &
     200                            + 6.0_wp * dist2(k,j,i)                                                &
    213201                            ) / 12.0_wp
    214202          ENDDO
    … …  
    219207
    220208!
    221 !-- Remove perturbations below topography (including one gridpoint above it
    222 !-- in order to allow for larger timesteps at the beginning of the simulation
    223 !-- (diffusion criterion))
     209!-- Remove perturbations below topography (including one gridpoint above it in order to allow for
     210!-- larger timesteps at the beginning of the simulation (diffusion criterion))
    224211    IF ( TRIM( topography ) /= 'flat' )  THEN
    225212       DO  i = nxlg, nxrg
    226213          DO  j = nysg, nyng
    227214             DO  k = nzb, nzb_max
    228                 dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                    &
    229                                   BTEST( wall_flags_total_0(k,j,i), flag_nr )  &
     215                dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                                        &
     216                                  BTEST( wall_flags_total_0(k,j,i), flag_nr )                      &
    230217                                    )
    231218             ENDDO
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |