Changeset 4488 for palm/trunk/SOURCE


Ignore:
Timestamp:
Apr 3, 2020 11:34:29 AM (4 years ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
8 edited

Legend:

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

    r4429 r4488  
    11!> @file advec_s_bc.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! 4429 2020-02-27 15:24:30Z raasch
    2729! bugfix: cpp-directives added for serial mode
    28 ! 
     30!
    2931! 4360 2020-01-07 11:25:50Z suehring
    3032! Corrected "Former revisions" section
    31 ! 
     33!
    3234! 3761 2019-02-25 15:31:42Z raasch
    3335! unused variables removed
    34 ! 
     36!
    3537! 3655 2019-01-07 16:51:22Z knoop
    3638! nopointer option removed
     
    4547!> Computation in individual steps for each of the three dimensions.
    4648!> Limiting assumptions:
    47 !> So far the scheme has been assuming equidistant grid spacing. As this is not
    48 !> the case in the stretched portion of the z-direction, there dzw(k) is used as
    49 !> a substitute for a constant grid length. This certainly causes incorrect
    50 !> results; however, it is hoped that they are not too apparent for weakly
    51 !> stretched grids.
     49!> So far the scheme has been assuming equidistant grid spacing. As this is not the case in the
     50!> stretched portion of the z-direction, there dzw(k) is used as a substitute for a constant grid
     51!> length. This certainly causes incorrect results; however, it is hoped that they are not too
     52!> apparent for weakly stretched grids.
    5253!> NOTE: This is a provisional, non-optimised version!
    53 !------------------------------------------------------------------------------!
    54 MODULE advec_s_bc_mod
    55  
     54!--------------------------------------------------------------------------------------------------!
     55 MODULE advec_s_bc_mod
     56
    5657
    5758    PRIVATE
     
    6465 CONTAINS
    6566
    66 !------------------------------------------------------------------------------!
     67!--------------------------------------------------------------------------------------------------!
    6768! Description:
    6869! ------------
    6970!> @todo Missing subroutine description.
    70 !------------------------------------------------------------------------------!
     71!--------------------------------------------------------------------------------------------------!
    7172    SUBROUTINE advec_s_bc( sk, sk_char )
    7273
    73        USE advection,                                                             &
     74       USE advection,                                                                              &
    7475           ONLY:  aex, bex, dex, eex
    7576
    76        USE arrays_3d,                                                             &
     77       USE arrays_3d,                                                                              &
    7778           ONLY:  d, ddzw, dzu, dzw, tend, u, v, w
    7879
    79        USE control_parameters,                                                    &
    80            ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, bc_s_t_val, ibc_pt_b, ibc_pt_t, &
    81                   ibc_q_t, ibc_s_t, message_string, pt_slope_offset,              &
    82                   sloping_surface, u_gtrans, v_gtrans
    83 
    84        USE cpulog,                                                                &
     80       USE control_parameters,                                                                     &
     81           ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, dt_3d, ibc_pt_b, ibc_pt_t, ibc_q_t, ibc_s_t,&
     82                  message_string, pt_slope_offset, sloping_surface, u_gtrans, v_gtrans
     83
     84       USE cpulog,                                                                                 &
    8585           ONLY:  cpu_log, log_point_s
    8686
    87        USE grid_variables,                                                        &
     87       USE grid_variables,                                                                         &
    8888           ONLY:  ddx, ddy
    8989
    90        USE indices,                                                               &
     90       USE indices,                                                                                &
    9191           ONLY:  nx, nxl, nxr, nyn, nys, nzb, nzt
    9292
     
    9595       USE pegrid
    9696
    97        USE statistics,                                                            &
     97       USE statistics,                                                                             &
    9898           ONLY:  rmask, statistic_regions, sums_wsts_bc_l
    9999
     
    112112       INTEGER(iwp) ::  type_xz_2 !<
    113113#endif
    114 
    115        REAL(wp) ::  cim    !<
    116        REAL(wp) ::  cimf   !<
    117        REAL(wp) ::  cip    !<
    118        REAL(wp) ::  cipf   !<
    119        REAL(wp) ::  d_new  !<
    120        REAL(wp) ::  denomi !< denominator
    121        REAL(wp) ::  fminus !<
    122        REAL(wp) ::  fplus  !<
    123        REAL(wp) ::  f2     !<
    124        REAL(wp) ::  f4     !<
    125        REAL(wp) ::  f8     !<
    126        REAL(wp) ::  f12    !<
    127        REAL(wp) ::  f24    !<
    128        REAL(wp) ::  f48    !<
    129        REAL(wp) ::  f1920  !<
    130        REAL(wp) ::  im     !<
    131        REAL(wp) ::  ip     !<
    132        REAL(wp) ::  m1n    !<
    133        REAL(wp) ::  m1z    !<
    134        REAL(wp) ::  m2     !<
    135        REAL(wp) ::  m3     !<
    136        REAL(wp) ::  numera !< numerator
    137        REAL(wp) ::  snenn  !<
    138        REAL(wp) ::  sterm  !<
    139        REAL(wp) ::  tendcy !<
    140        REAL(wp) ::  t1     !<
    141        REAL(wp) ::  t2     !<
     114       REAL(wp) ::  cim       !<
     115       REAL(wp) ::  cimf      !<
     116       REAL(wp) ::  cip       !<
     117       REAL(wp) ::  cipf      !<
     118       REAL(wp) ::  d_new     !<
     119       REAL(wp) ::  denomi    !< denominator
     120       REAL(wp) ::  fminus    !<
     121       REAL(wp) ::  fplus     !<
     122       REAL(wp) ::  f2        !<
     123       REAL(wp) ::  f4        !<
     124       REAL(wp) ::  f8        !<
     125       REAL(wp) ::  f12       !<
     126       REAL(wp) ::  f24       !<
     127       REAL(wp) ::  f48       !<
     128       REAL(wp) ::  f1920     !<
     129       REAL(wp) ::  im        !<
     130       REAL(wp) ::  ip        !<
     131       REAL(wp) ::  m1n       !<
     132       REAL(wp) ::  m1z       !<
     133       REAL(wp) ::  m2        !<
     134       REAL(wp) ::  m3        !<
     135       REAL(wp) ::  numera    !< numerator
     136       REAL(wp) ::  snenn     !<
     137       REAL(wp) ::  sterm     !<
     138       REAL(wp) ::  tendcy    !<
     139       REAL(wp) ::  t1        !<
     140       REAL(wp) ::  t2        !<
    142141
    143142       REAL(wp) ::  fmax(2)   !<
    144143       REAL(wp) ::  fmax_l(2) !<
    145        
    146        REAL(wp), DIMENSION(:,:,:), POINTER ::  sk
    147144
    148145       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a0   !<
     
    161158       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m1   !<
    162159       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sw   !<
    163        
     160
    164161       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sk_p !<
     162
     163       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !<
     164
    165165
    166166!
     
    197197!
    198198!--    Send left boundary, receive right boundary
    199        CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft,  0,      &
    200                           sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0,      &
     199       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft,  0,                       &
     200                          sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0,                       &
    201201                          comm2d, status, ierr )
    202202!
    203203!--    Send right boundary, receive left boundary
    204        CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1,      &
    205                           sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft,  1,      &
     204       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1,                       &
     205                          sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft,  1,                       &
    206206                          comm2d, status, ierr )
    207207       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
     
    217217
    218218!
    219 !--    In case of a sloping surface, the additional gridpoints in x-direction
    220 !--    of the temperature field at the left and right boundary of the total
    221 !--    domain must be adjusted by the temperature difference between this distance
     219!--    In case of a sloping surface, the additional gridpoints in x-direction of the temperature
     220!--    field at the left and right boundary of the total domain must be adjusted by the temperature
     221!--    difference between this distance
    222222       IF ( sloping_surface  .AND.  sk_char == 'pt' )  THEN
    223223          IF ( nxl ==  0 )  THEN
     
    241241          DO  j = nys, nyn
    242242             DO  k = nzb+1, nzt
    243                 numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
     243                numera  = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
    244244                denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    245245                fmax_l(1) = MAX( fmax_l(1) , numera )
     
    259259!
    260260!--    Allocate temporary arrays
    261        ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),          &
    262                  a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),         &
    263                  a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1),        &
    264                  imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1),        &
    265                  impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1),        &
    266                  ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1),        &
    267                  ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),          &
    268                  sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
     261       ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1), a2(nzb+1:nzt,nxl-1:nxr+1),&
     262                 a12(nzb+1:nzt,nxl-1:nxr+1),  a22(nzb+1:nzt,nxl-1:nxr+1),                          &
     263                 immb(nzb+1:nzt,nxl-1:nxr+1), imme(nzb+1:nzt,nxl-1:nxr+1),                         &
     264                 impb(nzb+1:nzt,nxl-1:nxr+1), impe(nzb+1:nzt,nxl-1:nxr+1),                         &
     265                 ipmb(nzb+1:nzt,nxl-1:nxr+1), ipme(nzb+1:nzt,nxl-1:nxr+1),                         &
     266                 ippb(nzb+1:nzt,nxl-1:nxr+1), ippe(nzb+1:nzt,nxl-1:nxr+1),                         &
     267                 m1(nzb+1:nzt,nxl-2:nxr+2),   sw(nzb+1:nzt,nxl-1:nxr+1)                            &
    269268               )
    270269       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    271270
    272271!
    273 !--    Initialise point of time measuring of the exponential portion (this would
    274 !--    not work if done locally within the loop)
     272!--    Initialise point of time measuring of the exponential portion (this would not work if done
     273!--    locally within the loop)
    275274       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
    276275       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
    277276
    278 !
    279 !--    Outer loop of all j
    280277       DO  j = nys, nyn
    281278
     
    285282             DO  k = nzb+1, nzt
    286283                a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
    287                 a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
    288                                                     + sk_p(k,j,i-1) )
    289                 a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
    290                             + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
    291                             + 9.0_wp * sk_p(k,j,i-2) ) * f1920
    292                 a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
    293                             - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
     284                a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
     285                a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)                   &
     286                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)                   &
     287                            + 9.0_wp * sk_p(k,j,i-2)                                               &
     288                          ) * f1920
     289                a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)                    &
     290                            - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)                     &
    294291                          ) * f48
    295                 a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
    296                             - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
    297                             - 3.0_wp * sk_p(k,j,i-2) ) * f48
     292                a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)                      &
     293                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)                      &
     294                            - 3.0_wp * sk_p(k,j,i-2)                                               &
     295                          ) * f48
    298296             ENDDO
    299297          ENDDO
     
    301299!
    302300!--       Fluxes using the Bott scheme
    303 !--       *VOCL LOOP,UNROLL(2)
    304301          DO  i = nxl, nxr
    305302             DO  k = nzb+1, nzt
     
    308305                cipf = 1.0_wp - 2.0_wp * cip
    309306                cimf = 1.0_wp - 2.0_wp * cim
    310                 ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
    311                        + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     307                ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                                       &
     308                       + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                                  &
    312309                       + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    313                 im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
    314                        - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     310                im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                                       &
     311                       - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                                  &
    315312                       + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    316313                ip   = MAX( ip, 0.0_wp )
     
    323320                cipf = 1.0_wp - 2.0_wp * cip
    324321                cimf = 1.0_wp - 2.0_wp * cim
    325                 ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
    326                        + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     322                ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                                       &
     323                       + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                                  &
    327324                       + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    328                 im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
    329                        - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     325                im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                                       &
     326                       - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                                  &
    330327                       + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    331328                ip   = MAX( ip, 0.0_wp )
     
    358355          DO  i = nxl-1, nxr+1
    359356             DO  k = nzb+1, nzt
    360                 m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
    361                      MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
     357                m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) / MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
    362358                IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
    363359
    364                 m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
    365                      MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
     360                m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) / MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
    366361                IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
    367362
     
    369364                t2 = 0.35_wp
    370365                IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
    371 
    372 !--             *VOCL STMT,IF(10)
    373                 IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
    374                      .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    375                      ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
    376                        m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
     366                IF ( m1(k,i-1) == 1.0_wp  .OR.  m1(k,i) == 1.0_wp  .OR.  m1(k,i+1) == 1.0_wp  .OR. &
     367                     m2 > t2  .OR.  m3 > t2  .OR.  ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0_wp .AND. &
     368                     m1(k,i) /= -1.0_wp .AND. m1(k,i+1) /= -1.0_wp )                               &
    377369                   )  sw(k,i) = 1.0_wp
    378370             ENDDO
     
    385377             DO  k = nzb+1, nzt
    386378
    387 !--             *VOCL STMT,IF(10)
    388379                IF ( sw(k,i) == 1.0_wp )  THEN
    389380                   snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
     
    397388                   cip =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    398389
    399                    ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
    400                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    401                                eex(ix) -                                          &
    402                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    403                                                                    )              &
     390                   ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                                     &
     391                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     392                               eex(ix) -                                                           &
     393                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     394                                                                   )                               &
    404395                                                             )
    405396                   IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
     
    416407                   cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    417408
    418                    imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
    419                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    420                                eex(ix) -                                          &
    421                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    422                                                                    )              &
     409                   imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                                     &
     410                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     411                               eex(ix) -                                                           &
     412                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     413                                                                   )                               &
    423414                                                             )
    424415                   IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
     
    426417                ENDIF
    427418
    428 !--             *VOCL STMT,IF(10)
    429419                IF ( sw(k,i+1) == 1.0_wp )  THEN
    430420                   snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
     
    438428                   cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
    439429
    440                    impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
    441                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    442                                eex(ix) -                                          &
    443                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    444                                                                    )              &
     430                   impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                                     &
     431                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     432                               eex(ix) -                                                           &
     433                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     434                                                                   )                               &
    445435                                                             )
    446436                   IF ( sterm == 0.0001_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
     
    448438                ENDIF
    449439
    450 !--             *VOCL STMT,IF(10)
    451440                IF ( sw(k,i-1) == 1.0_wp )  THEN
    452441                   snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
     
    460449                   cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
    461450
    462                    ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
    463                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    464                                eex(ix) -                                          &
    465                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    466                                                                    )              &
     451                   ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                                     &
     452                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     453                               eex(ix) -                                                           &
     454                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     455                                                                   )                               &
    467456                                                             )
    468457                   IF ( sterm == 0.0001_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
     
    478467          DO  i = nxl, nxr
    479468             DO  k = nzb+1, nzt
    480                 fplus  = ( 1.0_wp - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i)  &
    481                        - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
    482                 fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i)  &
    483                        - ( 1.0_wp - sw(k,i)   ) * immb(k,i) - sw(k,i)   * imme(k,i)
     469                fplus  = ( 1.0_wp - sw(k,i) )     * ippb(k,i) + sw(k,i)   * ippe(k,i)              &
     470                         - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
     471                fminus = ( 1.0_wp - sw(k,i-1) )   * ipmb(k,i) + sw(k,i-1) * ipme(k,i)              &
     472                         - ( 1.0_wp - sw(k,i) )  * immb(k,i) - sw(k,i)   * imme(k,i)
    484473                tendcy = fplus - fminus
    485474!
     
    490479!--             Density correction because of possible remaining divergences
    491480                d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
    492                 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
    493                               ( 1.0_wp + d_new )
     481                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / ( 1.0_wp + d_new )
    494482                d(k,j,i)  = d_new
    495483             ENDDO
     
    500488!
    501489!--    Deallocate temporary arrays
    502        DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    503                    ippb, ippe, m1, sw )
     490       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, ippb, ippe, m1, sw )
    504491
    505492
     
    508495#if defined( __parallel )
    509496       ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
    510        CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL,             &
     497       CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL,                              &
    511498                             type_xz_2, ierr )
    512499       CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
     
    514501!--    Send front boundary, receive rear boundary
    515502       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
    516        CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3),   1, type_xz_2, psouth, 0,       &
    517                           sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0,       &
     503       CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3),   1, type_xz_2, psouth, 0,                        &
     504                          sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0,                        &
    518505                          comm2d, status, ierr )
    519506!
    520507!--    Send rear boundary, receive front boundary
    521        CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1,       &
    522                           sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1,       &
     508       CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1,                        &
     509                          sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1,                        &
    523510                          comm2d, status, ierr )
    524511       CALL MPI_TYPE_FREE( type_xz_2, ierr )
     
    561548!
    562549!--    Allocate temporary arrays
    563        ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),          &
    564                  a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),         &
    565                  a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1),        &
    566                  imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1),        &
    567                  impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1),        &
    568                  ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1),        &
    569                  ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),          &
    570                  sw(nzb+1:nzt,nys-1:nyn+1)                                        &
     550       ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1), a2(nzb+1:nzt,nys-1:nyn+1),&
     551                 a12(nzb+1:nzt,nys-1:nyn+1),  a22(nzb+1:nzt,nys-1:nyn+1),                          &
     552                 immb(nzb+1:nzt,nys-1:nyn+1), imme(nzb+1:nzt,nys-1:nyn+1),                         &
     553                 impb(nzb+1:nzt,nys-1:nyn+1), impe(nzb+1:nzt,nys-1:nyn+1),                         &
     554                 ipmb(nzb+1:nzt,nys-1:nyn+1), ipme(nzb+1:nzt,nys-1:nyn+1),                         &
     555                 ippb(nzb+1:nzt,nys-1:nyn+1), ippe(nzb+1:nzt,nys-1:nyn+1),                         &
     556                 m1(nzb+1:nzt,nys-2:nyn+2),   sw(nzb+1:nzt,nys-1:nyn+1)                            &
    571557               )
    572558       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    573559
    574 !
    575 !--    Outer loop of all i
    576560       DO  i = nxl, nxr
    577561
     
    581565             DO  k = nzb+1, nzt
    582566                a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
    583                 a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
    584                                                     + sk_p(k,j-1,i) )
    585                 a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
    586                             + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
    587                             + 9.0_wp * sk_p(k,j-2,i) ) * f1920
    588                 a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
    589                             - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
     567                a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
     568                a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)                   &
     569                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)                   &
     570                            + 9.0_wp * sk_p(k,j-2,i)                                               &
     571                          ) * f1920
     572                a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)                    &
     573                            - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)                    &
    590574                          ) * f48
    591                 a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
    592                             - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
    593                             - 3.0_wp * sk_p(k,j-2,i) ) * f48
     575                a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)                      &
     576                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)                      &
     577                            - 3.0_wp * sk_p(k,j-2,i)                                               &
     578                          ) * f48
    594579             ENDDO
    595580          ENDDO
     
    597582!
    598583!--       Fluxes using the Bott scheme
    599 !--       *VOCL LOOP,UNROLL(2)
    600584          DO  j = nys, nyn
    601585             DO  k = nzb+1, nzt
     
    604588                cipf = 1.0_wp - 2.0_wp * cip
    605589                cimf = 1.0_wp - 2.0_wp * cim
    606                 ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
    607                        + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     590                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                                       &
     591                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                                  &
    608592                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    609                 im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
    610                        - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
     593                im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                                       &
     594                       - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                                  &
    611595                       + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    612596                ip   = MAX( ip, 0.0_wp )
     
    619603                cipf = 1.0_wp - 2.0_wp * cip
    620604                cimf = 1.0_wp - 2.0_wp * cim
    621                 ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
    622                        + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
     605                ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                                       &
     606                       + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                                  &
    623607                       + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    624                 im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
    625                        - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     608                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                                       &
     609                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                                  &
    626610                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    627611                ip   = MAX( ip, 0.0_wp )
     
    654638          DO  j = nys-1, nyn+1
    655639             DO  k = nzb+1, nzt
    656                 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
     640                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                                          &
    657641                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    658642                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
    659643
    660                 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
     644                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                                          &
    661645                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    662646                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     
    666650                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    667651
    668 !--             *VOCL STMT,IF(10)
    669                 IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
    670                      .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    671                      ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
    672                        m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
     652                IF ( m1(k,j-1) == 1.0_wp  .OR.  m1(k,j) == 1.0_wp  .OR.  m1(k,j+1) == 1.0_wp  .OR. &
     653                     m2 > t2  .OR.  m3 > t2  .OR.  ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0_wp .AND. &
     654                     m1(k,j) /= -1.0_wp .AND. m1(k,j+1) /= -1.0_wp )                               &
    673655                   )  sw(k,j) = 1.0_wp
    674656             ENDDO
     
    681663             DO  k = nzb+1, nzt
    682664
    683 !--             *VOCL STMT,IF(10)
    684665                IF ( sw(k,j) == 1.0_wp )  THEN
    685666                   snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
     
    693674                   cip =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    694675
    695                    ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
    696                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    697                                eex(ix) -                                          &
    698                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    699                                                                    )              &
     676                   ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                                     &
     677                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     678                               eex(ix) -                                                           &
     679                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     680                                                                   )                               &
    700681                                                             )
    701682                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     
    712693                   cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    713694
    714                    imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
    715                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    716                                eex(ix) -                                          &
    717                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    718                                                                    )              &
     695                   imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                                     &
     696                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     697                               eex(ix) -                                                           &
     698                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     699                                                                   )                               &
    719700                                                             )
    720701                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
     
    722703                ENDIF
    723704
    724 !--             *VOCL STMT,IF(10)
    725705                IF ( sw(k,j+1) == 1.0_wp )  THEN
    726706                   snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
     
    734714                   cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
    735715
    736                    impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
    737                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    738                                eex(ix) -                                          &
    739                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    740                                                                    )              &
     716                   impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                                     &
     717                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     718                               eex(ix) -                                                           &
     719                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     720                                                                   )                               &
    741721                                                             )
    742722                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
     
    744724                ENDIF
    745725
    746 !--             *VOCL STMT,IF(10)
    747726                IF ( sw(k,j-1) == 1.0_wp )  THEN
    748727                   snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
     
    756735                   cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
    757736
    758                    ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
    759                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    760                                eex(ix) -                                          &
    761                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    762                                                                    )              &
     737                   ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                                     &
     738                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     739                               eex(ix) -                                                           &
     740                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     741                                                                   )                               &
    763742                                                             )
    764743                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
     
    774753          DO  j = nys, nyn
    775754             DO  k = nzb+1, nzt
    776                 fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    777                        - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
    778                 fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
    779                        - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     755                fplus  = ( 1.0_wp - sw(k,j) )     * ippb(k,j) + sw(k,j)   * ippe(k,j)              &
     756                         - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
     757                fminus = ( 1.0_wp - sw(k,j-1) )   * ipmb(k,j) + sw(k,j-1) * ipme(k,j)              &
     758                         - ( 1.0_wp - sw(k,j) )  * immb(k,j) - sw(k,j)   * imme(k,j)
    780759                tendcy = fplus - fminus
    781760!
     
    786765!--             Density correction because of possible remaining divergences
    787766                d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
    788                 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
    789                               ( 1.0_wp + d_new )
    790                 d(k,j,i)  = d_new
     767                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / ( 1.0_wp + d_new )
     768                d(k,j,i) = d_new
    791769             ENDDO
    792770          ENDDO
    793771
    794772       ENDDO   ! End of the advection in y-direction
     773
     774
    795775       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
    796776       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
     
    798778!
    799779!--    Deallocate temporary arrays
    800        DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    801                    ippb, ippe, m1, sw )
     780       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, ippb, ippe, m1, sw )
    802781
    803782
     
    959938
    960939!
    961 !--       Cloud water content boundary condition at the bottom boundary:
    962 !--       Dirichlet (fixed surface rain water content).
     940!--       Cloud water content boundary condition at the bottom boundary: Dirichlet (fixed surface
     941!--       rain water content).
    963942          DO  i = nxl, nxr
    964943             DO  j = nys, nyn
     
    981960
    982961!
    983 !--       Rain water content boundary condition at the bottom boundary:
    984 !--       Dirichlet (fixed surface rain water content).
     962!--       Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface
     963!--       rain water content).
    985964          DO  i = nxl, nxr
    986965             DO  j = nys, nyn
     
    1003982
    1004983!
    1005 !--       Cloud drop concentration boundary condition at the bottom boundary:
    1006 !--       Dirichlet (fixed surface cloud drop concentration).
     984!--       Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
     985!--       surface cloud drop concentration).
    1007986          DO  i = nxl, nxr
    1008987             DO  j = nys, nyn
     
    10251004
    10261005!
    1027 !--       Rain drop concentration boundary condition at the bottom boundary:
    1028 !--       Dirichlet (fixed surface rain drop concentration).
     1006!--       Rain drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
     1007!--       surface rain drop concentration).
    10291008          DO  i = nxl, nxr
    10301009             DO  j = nys, nyn
     
    10601039       ELSE
    10611040
    1062           WRITE( message_string, * ) 'no vertical boundary condi',                &
    1063                                    'tion for variable "', sk_char, '"'
     1041          WRITE( message_string, * ) 'no vertical boundary condition for variable "', sk_char, '"'
    10641042          CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 )
    1065          
     1043
    10661044       ENDIF
    10671045
     
    10901068!
    10911069!--    Allocate temporary arrays
    1092        ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),                  &
    1093                  a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),                 &
    1094                  a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn),                &
    1095                  imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn),                &
    1096                  impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn),                &
    1097                  ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn),                &
    1098                  ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn),                &
    1099                  sw(nzb:nzt+1,nys:nyn)                                            &
     1070       ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),   a2(nzb:nzt+1,nys:nyn),          &
     1071                 a12(nzb:nzt+1,nys:nyn),  a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn),        &
     1072                 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), impe(nzb+1:nzt,nys:nyn),        &
     1073                 ipmb(nzb+1:nzt,nys:nyn), ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn),        &
     1074                 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), sw(nzb:nzt+1,nys:nyn)           &
    11001075               )
    11011076       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
    11021077
    1103 !
    1104 !--    Outer loop of all i
    11051078       DO  i = nxl, nxr
    11061079
     
    11101083             DO  k = nzb, nzt+1
    11111084                a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
    1112                 a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
    1113                                                     + sk_p(k-1,j,i) )
    1114                 a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
    1115                             + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
    1116                             + 9.0_wp * sk_p(k-2,j,i) ) * f1920
    1117                 a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
    1118                             - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
     1085                a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
     1086                a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)                   &
     1087                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)                   &
     1088                            + 9.0_wp * sk_p(k-2,j,i)                                               &
     1089                          ) * f1920
     1090                a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)                    &
     1091                            - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)                    &
    11191092                          ) * f48
    1120                 a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
    1121                             - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
    1122                             - 3.0_wp * sk_p(k-2,j,i) ) * f48
     1093                a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)                      &
     1094                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)                      &
     1095                            - 3.0_wp * sk_p(k-2,j,i)                                               &
     1096                          ) * f48
    11231097             ENDDO
    11241098          ENDDO
     
    11261100!
    11271101!--       Fluxes using the Bott scheme
    1128 !--       *VOCL LOOP,UNROLL(2)
    11291102          DO  j = nys, nyn
    11301103             DO  k = nzb+1, nzt
     
    11331106                cipf = 1.0_wp - 2.0_wp * cip
    11341107                cimf = 1.0_wp - 2.0_wp * cim
    1135                 ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
    1136                        + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
     1108                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                                       &
     1109                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                                  &
    11371110                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
    1138                 im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
    1139                        - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
     1111                im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                                       &
     1112                       - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                                  &
    11401113                       + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
    11411114                ip   = MAX( ip, 0.0_wp )
     
    11481121                cipf = 1.0_wp - 2.0_wp * cip
    11491122                cimf = 1.0_wp - 2.0_wp * cim
    1150                 ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
    1151                        + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
     1123                ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                                       &
     1124                       + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                                  &
    11521125                       + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
    1153                 im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
    1154                        - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
     1126                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                                       &
     1127                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                                  &
    11551128                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
    11561129                ip   = MAX( ip, 0.0_wp )
     
    11831156          DO  j = nys, nyn
    11841157             DO  k = nzb, nzt+1
    1185                 m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
     1158                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                                          &
    11861159                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
    11871160                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
    11881161
    1189                 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
     1162                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                                          &
    11901163                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
    11911164                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
     
    11951168                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
    11961169
    1197 !--             *VOCL STMT,IF(10)
    1198                 IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
    1199                      .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
    1200                      ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
    1201                        m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
     1170                IF ( m1(k-1,j) == 1.0_wp  .OR.  m1(k,j) == 1.0_wp  .OR.  m1(k+1,j) == 1.0_wp  .OR. &
     1171                     m2 > t2  .OR.  m3 > t2  .OR.  ( m1(k,j) > t1 .AND.  m1(k-1,j) /= -1.0_wp .AND.&
     1172                     m1(k,j) /= -1.0_wp .AND. m1(k+1,j) /= -1.0_wp )                               &
    12021173                   )  sw(k,j) = 1.0_wp
    12031174             ENDDO
     
    12101181             DO  k = nzb+1, nzt
    12111182
    1212 !--             *VOCL STMT,IF(10)
    12131183                IF ( sw(k,j) == 1.0_wp )  THEN
    12141184                   snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
     
    12221192                   cip =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    12231193
    1224                    ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
    1225                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1226                                eex(ix) -                                          &
    1227                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    1228                                                                    )              &
     1194                   ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                                     &
     1195                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     1196                               eex(ix) -                                                           &
     1197                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     1198                                                                   )                               &
    12291199                                                             )
    12301200                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
     
    12411211                   cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    12421212
    1243                    imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
    1244                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1245                                eex(ix) -                                          &
    1246                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
    1247                                                                    )              &
     1213                   imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                                     &
     1214                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     1215                               eex(ix) -                                                           &
     1216                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     1217                                                                   )                               &
    12481218                                                             )
    12491219                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
     
    12511221                ENDIF
    12521222
    1253 !--             *VOCL STMT,IF(10)
    12541223                IF ( sw(k+1,j) == 1.0_wp )  THEN
    12551224                   snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
     
    12631232                   cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
    12641233
    1265                    impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
    1266                                aex(ix) * cim + bex(ix) / dex(ix) * (              &
    1267                                eex(ix) -                                          &
    1268                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
    1269                                                                    )              &
     1234                   impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                                     &
     1235                               aex(ix) * cim + bex(ix) / dex(ix) * (                               &
     1236                               eex(ix) -                                                           &
     1237                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )                   &
     1238                                                                   )                               &
    12701239                                                             )
    12711240                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
     
    12731242                ENDIF
    12741243
    1275 !--             *VOCL STMT,IF(10)
    12761244                IF ( sw(k-1,j) == 1.0_wp )  THEN
    12771245                   snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
     
    12851253                   cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
    12861254
    1287                    ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
    1288                                aex(ix) * cip + bex(ix) / dex(ix) * (              &
    1289                                eex(ix) -                                          &
    1290                                EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
    1291                                                                    )              &
     1255                   ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                                     &
     1256                               aex(ix) * cip + bex(ix) / dex(ix) * (                               &
     1257                               eex(ix) -                                                           &
     1258                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )                   &
     1259                                                                   )                               &
    12921260                                                             )
    12931261                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
     
    13031271          DO  j = nys, nyn
    13041272             DO  k = nzb+1, nzt
    1305                 fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
    1306                        - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
    1307                 fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
    1308                        - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
     1273                fplus  = ( 1.0_wp - sw(k,j) )     * ippb(k,j) + sw(k,j)   * ippe(k,j)              &
     1274                         - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
     1275                fminus = ( 1.0_wp - sw(k-1,j) )   * ipmb(k,j) + sw(k-1,j) * ipme(k,j)              &
     1276                         - ( 1.0_wp - sw(k,j) )  * immb(k,j) - sw(k,j)   * imme(k,j)
    13091277                tendcy = fplus - fminus
    13101278!
     
    13151283!--             Density correction because of possible remaining divergences
    13161284                d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
    1317                 sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
    1318                               ( 1.0_wp + d_new )
     1285                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / ( 1.0_wp + d_new )
    13191286!
    13201287!--             Store heat flux for subsequent statistics output.
    1321 !--             array m1 is here used as temporary storage
     1288!--             Array m1 is here used as temporary storage
    13221289                m1(k,j) = fplus / dt_3d * dzw(k)
    13231290             ENDDO
     
    13251292
    13261293!
    1327 !--       Sum up heat flux in order to order to obtain horizontal averages
     1294!--       Sum up heat flux in order to obtain horizontal averages
    13281295          IF ( sk_char == 'pt' )  THEN
    13291296             DO  sr = 0, statistic_regions
    13301297                DO  j = nys, nyn
    13311298                   DO  k = nzb+1, nzt
    1332                       sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) +               &
    1333                                              m1(k,j) * rmask(j,i,sr)
     1299                      sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + m1(k,j) * rmask(j,i,sr)
    13341300                   ENDDO
    13351301                ENDDO
     
    13371303          ENDIF
    13381304
    1339        ENDDO   ! End of the advection in z-direction
     1305       ENDDO
     1306!
     1307!--    End of the advection in z-direction
     1308
    13401309       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
    13411310       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
     
    13431312!
    13441313!--    Deallocate temporary arrays
    1345        DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
    1346                    ippb, ippe, m1, sw )
     1314       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, ippb, ippe, m1, sw )
    13471315
    13481316!
  • palm/trunk/SOURCE/advec_s_pw.f90

    r4360 r4488  
    11!> @file advec_s_pw.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! 3927 2019-04-23 13:24:29Z raasch
    3032! pointer attribute removed from scalar 3d-array for performance reasons
    31 ! 
     33!
    3234! 3665 2019-01-10 08:28:24Z raasch
    3335! unused variables removed
    34 ! 
     36!
    3537! 3655 2019-01-07 16:51:22Z knoop
    3638! nopointer option removed
     
    4244! Description:
    4345! ------------
    44 !> Advection term for scalar variables using the Piacsek and Williams scheme
    45 !> (form C3). Contrary to PW itself, for reasons of accuracy their scheme is
    46 !> slightly modified as follows: the values of those scalars that are used for
    47 !> the computation of the flux divergence are reduced by the value of the
    48 !> relevant scalar at the location where the difference is computed (sk(k,j,i)).
     46!> Advection term for scalar variables using the Piacsek and Williams scheme (form C3). Contrary to
     47!> PW itself, for reasons of accuracy their scheme is slightly modified as follows: the values of
     48!> those scalars that are used for the computation of the flux divergence are reduced by the value
     49!> of the relevant scalar at the location where the difference is computed (sk(k,j,i)).
    4950!> NOTE: at the first grid point above the surface computation still takes place!
    50 !------------------------------------------------------------------------------!
     51!--------------------------------------------------------------------------------------------------!
    5152 MODULE advec_s_pw_mod
    52  
     53
    5354
    5455    PRIVATE
     
    5960       MODULE PROCEDURE advec_s_pw_ij
    6061    END INTERFACE
    61  
     62
    6263 CONTAINS
    6364
    6465
    65 !------------------------------------------------------------------------------!
     66!--------------------------------------------------------------------------------------------------!
    6667! Description:
    6768! ------------
    6869!> Call for all grid points
    69 !------------------------------------------------------------------------------!
    70     SUBROUTINE advec_s_pw( sk )
     70!--------------------------------------------------------------------------------------------------!
     71 SUBROUTINE advec_s_pw( sk )
    7172
    72        USE arrays_3d,                                                          &
    73            ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
     73    USE arrays_3d,                                                                                 &
     74        ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
    7475
    75        USE control_parameters,                                                 &
    76            ONLY:  u_gtrans, v_gtrans
     76    USE control_parameters,                                                                        &
     77        ONLY:  u_gtrans, v_gtrans
    7778
    78        USE grid_variables,                                                     &
    79            ONLY:  ddx, ddy
     79    USE grid_variables,                                                                            &
     80        ONLY:  ddx, ddy
    8081
    81        USE indices,                                                            &
    82            ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
     82    USE indices,                                                                                   &
     83        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
    8384
    84        USE kinds
     85    USE kinds
    8586
    8687
    87        IMPLICIT NONE
     88    IMPLICIT NONE
    8889
    89        INTEGER(iwp) ::  i !< grid index along x-direction
    90        INTEGER(iwp) ::  j !< grid index along y-direction
    91        INTEGER(iwp) ::  k !< grid index along z-direction
     90    INTEGER(iwp) ::  i !< grid index along x-direction
     91    INTEGER(iwp) ::  j !< grid index along y-direction
     92    INTEGER(iwp) ::  k !< grid index along z-direction
    9293
    93        REAL(wp)     ::  gu  !< local additional advective velocity
    94        REAL(wp)     ::  gv  !< local additional advective velocity
     94    REAL(wp)     ::  gu  !< local additional advective velocity
     95    REAL(wp)     ::  gv  !< local additional advective velocity
    9596
    96        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
    97  
     97    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
    9898
    99        DO  i = nxl, nxr
    100           DO  j = nys, nyn
    101              DO  k = nzb+1, nzt
     99
     100    DO  i = nxl, nxr
     101       DO  j = nys, nyn
     102          DO  k = nzb+1, nzt
    102103
    103104!
    104 !--             Galilean transformation + Stokes drift velocity (ocean only)
    105                 gu = u_gtrans - u_stokes_zu(k)
    106                 gv = v_gtrans - v_stokes_zu(k)
     105!--          Galilean transformation + Stokes drift velocity (ocean only)
     106             gu = u_gtrans - u_stokes_zu(k)
     107             gv = v_gtrans - v_stokes_zu(k)
    107108
    108                 tend(k,j,i) = tend(k,j,i) +                                    &
    109                                      ( -0.5_wp * ( ( u(k,j,i+1) - gu       ) * &
    110                                                    ( sk(k,j,i+1) - sk(k,j,i) ) &
    111                                                  - ( u(k,j,i)   - gu       ) * &
    112                                                    ( sk(k,j,i-1) - sk(k,j,i) ) &
    113                                                  ) * ddx                       &
    114                                        -0.5_wp * ( ( v(k,j+1,i) - gv       ) * &
    115                                                    ( sk(k,j+1,i) - sk(k,j,i) ) &
    116                                                  - ( v(k,j,i)   - gv       ) * &
    117                                                    ( sk(k,j-1,i) - sk(k,j,i) ) &
    118                                                  ) * ddy                       &
    119                                        -         (   w(k,j,i)                * &
    120                                                    ( sk(k+1,j,i) - sk(k,j,i) ) &
    121                                                  -   w(k-1,j,i)              * &
    122                                                    ( sk(k-1,j,i) - sk(k,j,i) ) &
    123                                                  ) * dd2zu(k)                  &
    124                                       )
    125              ENDDO
     109             tend(k,j,i) = tend(k,j,i) +                                                           &
     110                                  ( -0.5_wp * ( ( u(k,j,i+1)  - gu        ) *                      &
     111                                                ( sk(k,j,i+1) - sk(k,j,i) )                        &
     112                                              - ( u(k,j,i)    - gu        ) *                      &
     113                                                ( sk(k,j,i-1) - sk(k,j,i) )                        &
     114                                              ) * ddx                                              &
     115                                    -0.5_wp * ( ( v(k,j+1,i)  - gv        ) *                      &
     116                                                ( sk(k,j+1,i) - sk(k,j,i) )                        &
     117                                              - ( v(k,j,i)    - gv        ) *                      &
     118                                                ( sk(k,j-1,i) - sk(k,j,i) )                        &
     119                                              ) * ddy                                              &
     120                                    -         (   w(k,j,i)                  *                      &
     121                                                ( sk(k+1,j,i) - sk(k,j,i) )                        &
     122                                              -   w(k-1,j,i)                *                      &
     123                                                ( sk(k-1,j,i) - sk(k,j,i) )                        &
     124                                              ) * dd2zu(k)                                         &
     125                                   )
    126126          ENDDO
    127127       ENDDO
     128    ENDDO
    128129
    129     END SUBROUTINE advec_s_pw
     130 END SUBROUTINE advec_s_pw
    130131
    131132
    132 !------------------------------------------------------------------------------!
     133!--------------------------------------------------------------------------------------------------!
    133134! Description:
    134135! ------------
    135136!> Call for grid point i,j
    136 !------------------------------------------------------------------------------!
    137     SUBROUTINE advec_s_pw_ij( i, j, sk )
     137!--------------------------------------------------------------------------------------------------!
     138 SUBROUTINE advec_s_pw_ij( i, j, sk )
    138139
    139        USE arrays_3d,                                                          &
    140            ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
     140    USE arrays_3d,                                                                                 &
     141        ONLY:  dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w
    141142
    142        USE control_parameters,                                                 &
    143            ONLY:  u_gtrans, v_gtrans
     143    USE control_parameters,                                                                        &
     144        ONLY:  u_gtrans, v_gtrans
    144145
    145        USE grid_variables,                                                     &
    146            ONLY:  ddx, ddy
     146    USE grid_variables,                                                                            &
     147        ONLY:  ddx, ddy
    147148
    148        USE indices,                                                            &
    149            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
     149    USE indices,                                                                                   &
     150        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
    150151
    151        USE kinds
     152    USE kinds
    152153
    153154
    154        IMPLICIT NONE
     155    IMPLICIT NONE
    155156
    156        INTEGER(iwp) ::  i !< grid index along x-direction
    157        INTEGER(iwp) ::  j !< grid index along y-direction
    158        INTEGER(iwp) ::  k !< grid index along z-direction
     157    INTEGER(iwp) ::  i !< grid index along x-direction
     158    INTEGER(iwp) ::  j !< grid index along y-direction
     159    INTEGER(iwp) ::  k !< grid index along z-direction
    159160
    160        REAL(wp)     ::  gu  !< local additional advective velocity
    161        REAL(wp)     ::  gv  !< local additional advective velocity
     161    REAL(wp)     ::  gu  !< local additional advective velocity
     162    REAL(wp)     ::  gv  !< local additional advective velocity
    162163
    163        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
     164    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
    164165
    165166
    166        DO  k = nzb+1, nzt
     167    DO  k = nzb+1, nzt
    167168
    168169!
    169 !--       Galilean transformation + Stokes drift velocity (ocean only)
    170           gu = u_gtrans - u_stokes_zu(k)
    171           gv = v_gtrans - v_stokes_zu(k)
     170!--    Galilean transformation + Stokes drift velocity (ocean only)
     171       gu = u_gtrans - u_stokes_zu(k)
     172       gv = v_gtrans - v_stokes_zu(k)
    172173
    173           tend(k,j,i) = tend(k,j,i) +                                          &
    174                                     ( -0.5_wp * ( ( u(k,j,i+1) - gu        ) * &
    175                                                   ( sk(k,j,i+1) - sk(k,j,i) )  &
    176                                                 - ( u(k,j,i)   - gu        ) * &
    177                                                   ( sk(k,j,i-1) - sk(k,j,i) )  &
    178                                                 ) * ddx                        &
    179                                       -0.5_wp * ( ( v(k,j+1,i) - gv       ) *  &
    180                                                   ( sk(k,j+1,i) - sk(k,j,i) )  &
    181                                                 - ( v(k,j,i)   - gv       ) *  &
    182                                                   ( sk(k,j-1,i) - sk(k,j,i) )  &
    183                                                 ) * ddy                        &
    184                                       -         (   w(k,j,i)                *  &
    185                                                   ( sk(k+1,j,i) - sk(k,j,i) )  &
    186                                                 -   w(k-1,j,i)              *  &
    187                                                   ( sk(k-1,j,i) - sk(k,j,i) )  &
    188                                                 ) * dd2zu(k)                   &
    189                                     )
    190        ENDDO
     174       tend(k,j,i) = tend(k,j,i) +                                                                 &
     175                                 ( -0.5_wp * ( ( u(k,j,i+1)  - gu        ) *                      &
     176                                               ( sk(k,j,i+1) - sk(k,j,i) )                         &
     177                                             - ( u(k,j,i)    - gu        ) *                      &
     178                                               ( sk(k,j,i-1) - sk(k,j,i) )                         &
     179                                             ) * ddx                                               &
     180                                   -0.5_wp * ( ( v(k,j+1,i)  - gv        ) *                       &
     181                                               ( sk(k,j+1,i) - sk(k,j,i) )                         &
     182                                             - ( v(k,j,i)    - gv        ) *                       &
     183                                               ( sk(k,j-1,i) - sk(k,j,i) )                         &
     184                                             ) * ddy                                               &
     185                                   -         (   w(k,j,i)                  *                       &
     186                                               ( sk(k+1,j,i) - sk(k,j,i) )                         &
     187                                             -   w(k-1,j,i)                *                       &
     188                                               ( sk(k-1,j,i) - sk(k,j,i) )                         &
     189                                             ) * dd2zu(k)                                          &
     190                                 )
     191    ENDDO
    191192
    192     END SUBROUTINE advec_s_pw_ij
     193 END SUBROUTINE advec_s_pw_ij
    193194
    194195 END MODULE advec_s_pw_mod
  • palm/trunk/SOURCE/advec_s_up.f90

    r4360 r4488  
    11!> @file advec_s_up.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! 3927 2019-04-23 13:24:29Z raasch
    3032! pointer attribute removed from scalar 3d-array for performance reasons
    31 ! 
     33!
    3234! 3665 2019-01-10 08:28:24Z raasch
    3335! unused variables removed
    34 ! 
     36!
    3537! 3655 2019-01-07 16:51:22Z knoop
    3638! nopointer option removed
     
    4547!> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0!
    4648!>       The same problem occurs for all topography boundaries!
    47 !------------------------------------------------------------------------------!
     49!--------------------------------------------------------------------------------------------------!
    4850 MODULE advec_s_up_mod
    49  
     51
    5052
    5153    PRIVATE
     
    6062
    6163
    62 !------------------------------------------------------------------------------!
     64!--------------------------------------------------------------------------------------------------!
    6365! Description:
    6466! ------------
    6567!> Call for all grid points
    66 !------------------------------------------------------------------------------!
    67     SUBROUTINE advec_s_up( sk )
     68!--------------------------------------------------------------------------------------------------!
     69 SUBROUTINE advec_s_up( sk )
    6870
    69        USE arrays_3d,                                                          &
    70            ONLY:  ddzu, tend, u, v, w
     71    USE arrays_3d,                                                                                 &
     72        ONLY:  ddzu, tend, u, v, w
    7173
    72        USE control_parameters,                                                 &
    73            ONLY:  u_gtrans, v_gtrans
     74    USE control_parameters,                                                                        &
     75        ONLY:  u_gtrans, v_gtrans
    7476
    75        USE grid_variables,                                                     &
    76            ONLY:  ddx, ddy
     77    USE grid_variables,                                                                            &
     78        ONLY:  ddx, ddy
    7779
    78        USE indices,                                                            &
    79            ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
     80    USE indices,                                                                                   &
     81        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
    8082
    81        USE kinds
     83    USE kinds
    8284
    8385
    84        IMPLICIT NONE
     86    IMPLICIT NONE
    8587
    86        INTEGER(iwp) ::  i !< grid index along x-direction
    87        INTEGER(iwp) ::  j !< grid index along y-direction
    88        INTEGER(iwp) ::  k !< grid index along z-direction
     88    INTEGER(iwp) ::  i !< grid index along x-direction
     89    INTEGER(iwp) ::  j !< grid index along y-direction
     90    INTEGER(iwp) ::  k !< grid index along z-direction
    8991
    90        REAL(wp) ::  ukomp !< advection velocity along x-direction
    91        REAL(wp) ::  vkomp !< advection velocity along y-direction
    92        REAL(wp) ::  wkomp !< advection velocity along z-direction
     92    REAL(wp) ::  ukomp !< advection velocity along x-direction
     93    REAL(wp) ::  vkomp !< advection velocity along y-direction
     94    REAL(wp) ::  wkomp !< advection velocity along z-direction
    9395
    94        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
     96    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
    9597
    9698
    97        DO  i = nxl, nxr
    98           DO  j = nys, nyn
    99              DO  k = nzb+1, nzt
     99    DO  i = nxl, nxr
     100       DO  j = nys, nyn
     101          DO  k = nzb+1, nzt
    100102!
    101 !--             x-direction
    102                 ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
    103                 IF ( ukomp > 0.0_wp )  THEN
    104                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    105                                          ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
    106                 ELSE
    107                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    108                                          ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
    109                 ENDIF
     103!--          x-direction
     104             ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
     105             IF ( ukomp > 0.0_wp )  THEN
     106                tend(k,j,i) = tend(k,j,i) - ukomp * ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
     107             ELSE
     108                tend(k,j,i) = tend(k,j,i) - ukomp * ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
     109             ENDIF
    110110!
    111 !--             y-direction
    112                 vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
    113                 IF ( vkomp > 0.0_wp )  THEN
    114                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    115                                          ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
    116                 ELSE
    117                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    118                                          ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
    119                 ENDIF
     111!--          y-direction
     112             vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
     113             IF ( vkomp > 0.0_wp )  THEN
     114                tend(k,j,i) = tend(k,j,i) - vkomp * ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
     115             ELSE
     116                tend(k,j,i) = tend(k,j,i) - vkomp * ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
     117             ENDIF
    120118!
    121 !--             z-direction
    122                 wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
    123                 IF ( wkomp > 0.0_wp )  THEN
    124                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    125                                          ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
    126                 ELSE
    127                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    128                                          ( sk(k+1,j,i)-sk(k,j,i) ) * ddzu(k+1)
    129                 ENDIF
     119!--          z-direction
     120             wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     121             IF ( wkomp > 0.0_wp )  THEN
     122                tend(k,j,i) = tend(k,j,i) - wkomp * ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
     123             ELSE
     124                tend(k,j,i) = tend(k,j,i) - wkomp * ( sk(k+1,j,i) - sk(k,j,i) ) * ddzu(k+1)
     125             ENDIF
    130126
    131              ENDDO
    132127          ENDDO
    133128       ENDDO
     129    ENDDO
    134130
    135     END SUBROUTINE advec_s_up
     131 END SUBROUTINE advec_s_up
    136132
    137133
    138 !------------------------------------------------------------------------------!
     134!--------------------------------------------------------------------------------------------------!
    139135! Description:
    140136! ------------
    141137!> Call for grid point i,j
    142 !------------------------------------------------------------------------------!
    143     SUBROUTINE advec_s_up_ij( i, j, sk )
     138!--------------------------------------------------------------------------------------------------!
     139 SUBROUTINE advec_s_up_ij( i, j, sk )
    144140
    145        USE arrays_3d,                                                          &
    146            ONLY:  ddzu, tend, u, v, w
     141    USE arrays_3d,                                                                                 &
     142        ONLY:  ddzu, tend, u, v, w
    147143
    148        USE control_parameters,                                                 &
    149            ONLY:  u_gtrans, v_gtrans
     144    USE control_parameters,                                                                        &
     145        ONLY:  u_gtrans, v_gtrans
    150146
    151        USE grid_variables,                                                     &
    152            ONLY:  ddx, ddy
     147    USE grid_variables,                                                                            &
     148        ONLY:  ddx, ddy
    153149
    154        USE indices,                                                            &
    155            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
     150    USE indices,                                                                                   &
     151        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
    156152
    157        USE kinds
     153    USE kinds
    158154
    159155
    160        IMPLICIT NONE
     156    IMPLICIT NONE
    161157
    162        INTEGER(iwp) ::  i !< grid index along x-direction
    163        INTEGER(iwp) ::  j !< grid index along y-direction
    164        INTEGER(iwp) ::  k !< grid index along z-direction
     158    INTEGER(iwp) ::  i !< grid index along x-direction
     159    INTEGER(iwp) ::  j !< grid index along y-direction
     160    INTEGER(iwp) ::  k !< grid index along z-direction
    165161
    166        REAL(wp) ::  ukomp !< advection velocity along x-direction
    167        REAL(wp) ::  vkomp !< advection velocity along y-direction
    168        REAL(wp) ::  wkomp !< advection velocity along z-direction
    169        
    170        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
     162    REAL(wp) ::  ukomp !< advection velocity along x-direction
     163    REAL(wp) ::  vkomp !< advection velocity along y-direction
     164    REAL(wp) ::  wkomp !< advection velocity along z-direction
     165
     166    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
    171167
    172168
    173        DO  k = nzb+1, nzt
     169    DO  k = nzb+1, nzt
    174170!
    175 !--       x-direction
    176           ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
    177           IF ( ukomp > 0.0_wp )  THEN
    178              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    179                                          ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
    180           ELSE
    181              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    182                                          ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
    183           ENDIF
     171!--    x-direction
     172       ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
     173       IF ( ukomp > 0.0_wp )  THEN
     174          tend(k,j,i) = tend(k,j,i) - ukomp * ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
     175       ELSE
     176          tend(k,j,i) = tend(k,j,i) - ukomp * ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
     177       ENDIF
    184178!
    185 !--       y-direction
    186           vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
    187           IF ( vkomp > 0.0_wp )  THEN
    188              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    189                                          ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
    190           ELSE
    191              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    192                                          ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
    193           ENDIF
     179!--    y-direction
     180       vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
     181       IF ( vkomp > 0.0_wp )  THEN
     182          tend(k,j,i) = tend(k,j,i) - vkomp * ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
     183       ELSE
     184          tend(k,j,i) = tend(k,j,i) - vkomp * ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
     185       ENDIF
    194186!
    195 !--       z-direction
    196           wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
    197           IF ( wkomp > 0.0_wp )  THEN
    198              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    199                                          ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
    200           ELSE
    201              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    202                                          ( sk(k+1,j,i)-sk(k,j,i) ) * ddzu(k+1)
    203           ENDIF
     187!--    z-direction
     188       wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
     189       IF ( wkomp > 0.0_wp )  THEN
     190          tend(k,j,i) = tend(k,j,i) - wkomp * ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
     191       ELSE
     192          tend(k,j,i) = tend(k,j,i) - wkomp * ( sk(k+1,j,i) - sk(k,j,i) ) * ddzu(k+1)
     193       ENDIF
    204194
    205        ENDDO
     195    ENDDO
    206196
    207     END SUBROUTINE advec_s_up_ij
     197 END SUBROUTINE advec_s_up_ij
    208198
    209199 END MODULE advec_s_up_mod
  • palm/trunk/SOURCE/advec_u_pw.f90

    r4360 r4488  
    11!> @file advec_u_pw.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! variables documented
     
    3739! ------------
    3840!> Advection term for u velocity-component using Piacsek and Williams.
    39 !> Vertical advection at the first grid point above the surface is done with
    40 !> normal centred differences, because otherwise no information from the surface
    41 !> would be communicated upwards due to w=0 at K=nzb.
    42 !------------------------------------------------------------------------------!
     41!> Vertical advection at the first grid point above the surface is done with normal centred
     42!> differences, because otherwise no information from the surface would be communicated upwards due
     43!> to w=0 at K=nzb.
     44!--------------------------------------------------------------------------------------------------!
    4345 MODULE advec_u_pw_mod
    44  
     46
    4547
    4648    PRIVATE
     
    5153       MODULE PROCEDURE advec_u_pw_ij
    5254    END INTERFACE advec_u_pw
    53  
     55
    5456 CONTAINS
    5557
    5658
    57 !------------------------------------------------------------------------------!
     59!--------------------------------------------------------------------------------------------------!
    5860! Description:
    5961! ------------
    6062!> Call for all grid points
    61 !------------------------------------------------------------------------------!
    62     SUBROUTINE advec_u_pw
     63!--------------------------------------------------------------------------------------------------!
     64 SUBROUTINE advec_u_pw
    6365
    64        USE arrays_3d,                                                          &
    65            ONLY:  ddzw, tend, u, v, w
     66    USE arrays_3d,                                                                                 &
     67        ONLY:  ddzw, tend, u, v, w
    6668
    67        USE control_parameters,                                                 &
    68            ONLY:  u_gtrans, v_gtrans
     69    USE control_parameters,                                                                        &
     70        ONLY:  u_gtrans, v_gtrans
    6971
    70        USE grid_variables,                                                     &
    71            ONLY:  ddx, ddy
     72    USE grid_variables,                                                                            &
     73        ONLY:  ddx, ddy
    7274
    73        USE indices,                                                            &
    74            ONLY:  nxlu, nxr, nyn, nys, nzb, nzt
     75    USE indices,                                                                                   &
     76        ONLY:  nxlu, nxr, nyn, nys, nzb, nzt
    7577
    76        USE kinds
     78    USE kinds
    7779
    7880
    79        IMPLICIT NONE
     81    IMPLICIT NONE
    8082
    81        INTEGER(iwp) ::  i !< grid index along x-direction
    82        INTEGER(iwp) ::  j !< grid index along y-direction
    83        INTEGER(iwp) ::  k !< grid index along z-direction
    84        
    85        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    86        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
    87  
    88        gu = 2.0_wp * u_gtrans
    89        gv = 2.0_wp * v_gtrans
    90        DO  i = nxlu, nxr
    91           DO  j = nys, nyn
    92              DO  k = nzb+1, nzt
    93                 tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    94                          ( u(k,j,i+1) * ( u(k,j,i+1) + u(k,j,i) - gu )         &
    95                          - u(k,j,i-1) * ( u(k,j,i) + u(k,j,i-1) - gu ) ) * ddx &
    96                        + ( u(k,j+1,i) * ( v(k,j+1,i) + v(k,j+1,i-1) - gv )     &
    97                          - u(k,j-1,i) * ( v(k,j,i) + v(k,j,i-1) - gv ) ) * ddy &
    98                        + ( u(k+1,j,i) * ( w(k,j,i) + w(k,j,i-1) )              &
    99                          - u(k-1,j,i) * ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    100                                                                   * ddzw(k)    &
    101                                                       )
    102              ENDDO
     83    INTEGER(iwp) ::  i !< grid index along x-direction
     84    INTEGER(iwp) ::  j !< grid index along y-direction
     85    INTEGER(iwp) ::  k !< grid index along z-direction
     86
     87    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     88    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
     89
     90    gu = 2.0_wp * u_gtrans
     91    gv = 2.0_wp * v_gtrans
     92    DO  i = nxlu, nxr
     93       DO  j = nys, nyn
     94          DO  k = nzb+1, nzt
     95             tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                               &
     96                           ( u(k,j,i+1)   * ( u(k,j,i+1) + u(k,j,i) - gu )                         &
     97                           - u(k,j,i-1)   * ( u(k,j,i)   + u(k,j,i-1) - gu ) ) * ddx               &
     98                           + ( u(k,j+1,i) * ( v(k,j+1,i) + v(k,j+1,i-1) - gv )                     &
     99                           - u(k,j-1,i)   * ( v(k,j,i)   + v(k,j,i-1) - gv ) ) * ddy               &
     100                           + ( u(k+1,j,i) * ( w(k,j,i)   + w(k,j,i-1) )                            &
     101                           - u(k-1,j,i)   * ( w(k-1,j,i) + w(k-1,j,i-1) ) )    * ddzw(k)           &
     102                                                   )
    103103          ENDDO
    104104       ENDDO
     105    ENDDO
    105106
    106     END SUBROUTINE advec_u_pw
     107 END SUBROUTINE advec_u_pw
    107108
    108109
    109 !------------------------------------------------------------------------------!
     110!--------------------------------------------------------------------------------------------------!
    110111! Description:
    111112! ------------
    112113!> Call for grid point i,j
    113 !------------------------------------------------------------------------------!
    114     SUBROUTINE advec_u_pw_ij( i, j )
     114!--------------------------------------------------------------------------------------------------!
     115 SUBROUTINE advec_u_pw_ij( i, j )
    115116
    116        USE arrays_3d,                                                          &
    117            ONLY:  ddzw, tend, u, v, w
     117    USE arrays_3d,                                                                                 &
     118        ONLY:  ddzw, tend, u, v, w
    118119
    119        USE control_parameters,                                                 &
    120            ONLY:  u_gtrans, v_gtrans
     120    USE control_parameters,                                                                        &
     121        ONLY:  u_gtrans, v_gtrans
    121122
    122        USE grid_variables,                                                     &
    123            ONLY:  ddx, ddy
     123    USE grid_variables,                                                                            &
     124        ONLY:  ddx, ddy
    124125
    125        USE indices,                                                            &
    126            ONLY:  nzb, nzt
     126    USE indices,                                                                                   &
     127        ONLY:  nzb, nzt
    127128
    128        USE kinds
     129    USE kinds
    129130
    130131
    131        IMPLICIT NONE
     132    IMPLICIT NONE
    132133
    133        INTEGER(iwp) ::  i !< grid index along x-direction
    134        INTEGER(iwp) ::  j !< grid index along y-direction
    135        INTEGER(iwp) ::  k !< grid index along z-direction
    136        
    137        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    138        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
     134    INTEGER(iwp) ::  i !< grid index along x-direction
     135    INTEGER(iwp) ::  j !< grid index along y-direction
     136    INTEGER(iwp) ::  k !< grid index along z-direction
    139137
    140        gu = 2.0_wp * u_gtrans
    141        gv = 2.0_wp * v_gtrans
    142        DO  k = nzb+1, nzt
    143           tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    144                          ( u(k,j,i+1) * ( u(k,j,i+1) + u(k,j,i) - gu )         &
    145                          - u(k,j,i-1) * ( u(k,j,i) + u(k,j,i-1) - gu ) ) * ddx &
    146                        + ( u(k,j+1,i) * ( v(k,j+1,i) + v(k,j+1,i-1) - gv )     &
    147                          - u(k,j-1,i) * ( v(k,j,i) + v(k,j,i-1) - gv ) ) * ddy &
    148                        + ( u(k+1,j,i) * ( w(k,j,i) + w(k,j,i-1) )              &
    149                          - u(k-1,j,i) * ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    150                                                                   * ddzw(k)    &
    151                                                 )
    152        ENDDO
     138    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     139    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
    153140
    154     END SUBROUTINE advec_u_pw_ij
     141    gu = 2.0_wp * u_gtrans
     142    gv = 2.0_wp * v_gtrans
     143    DO  k = nzb+1, nzt
     144       tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                                     &
     145                     ( u(k,j,i+1)   * ( u(k,j,i+1) + u(k,j,i) - gu )                               &
     146                     - u(k,j,i-1)   * ( u(k,j,i)   + u(k,j,i-1) - gu ) ) * ddx                     &
     147                     + ( u(k,j+1,i) * ( v(k,j+1,i) + v(k,j+1,i-1) - gv )                           &
     148                     - u(k,j-1,i)   * ( v(k,j,i)   + v(k,j,i-1) - gv ) ) * ddy                     &
     149                     + ( u(k+1,j,i) * ( w(k,j,i)   + w(k,j,i-1) )                                  &
     150                     - u(k-1,j,i)   * ( w(k-1,j,i) + w(k-1,j,i-1) ) )    * ddzw(k)                 &
     151                                             )
     152    ENDDO
     153
     154 END SUBROUTINE advec_u_pw_ij
    155155
    156156 END MODULE advec_u_pw_mod
  • palm/trunk/SOURCE/advec_u_up.f90

    r4360 r4488  
    11!> @file advec_u_up.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! variables documented
     
    3941!> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0!
    4042!>       The same problem occurs for all topography boundaries!
    41 !------------------------------------------------------------------------------!
     43!--------------------------------------------------------------------------------------------------!
    4244 MODULE advec_u_up_mod
    43  
     45
    4446
    4547    PRIVATE
     
    5456
    5557
    56 !------------------------------------------------------------------------------!
     58!--------------------------------------------------------------------------------------------------!
    5759! Description:
    5860! ------------
    5961!> Call for all grid points
    60 !------------------------------------------------------------------------------!
    61     SUBROUTINE advec_u_up
     62!--------------------------------------------------------------------------------------------------!
     63 SUBROUTINE advec_u_up
    6264
    63        USE arrays_3d,                                                          &
    64            ONLY:  ddzu, tend, u, v, w
     65    USE arrays_3d,                                                                                 &
     66        ONLY:  ddzu, tend, u, v, w
    6567
    66        USE control_parameters,                                                 &
    67            ONLY:  u_gtrans, v_gtrans
     68    USE control_parameters,                                                                        &
     69        ONLY:  u_gtrans, v_gtrans
    6870
    69        USE grid_variables,                                                     &
    70            ONLY:  ddx, ddy
     71    USE grid_variables,                                                                            &
     72        ONLY:  ddx, ddy
    7173
    72        USE indices,                                                            &
    73            ONLY:  nxlu, nxr, nyn, nys, nzb, nzt
     74    USE indices,                                                                                   &
     75        ONLY:  nxlu, nxr, nyn, nys, nzb, nzt
    7476
    75        USE kinds
     77    USE kinds
    7678
    7779
    78        IMPLICIT NONE
     80    IMPLICIT NONE
    7981
    80        INTEGER(iwp) ::  i !< grid index along x-direction
    81        INTEGER(iwp) ::  j !< grid index along y-direction
    82        INTEGER(iwp) ::  k !< grid index along z-direction
     82    INTEGER(iwp) ::  i !< grid index along x-direction
     83    INTEGER(iwp) ::  j !< grid index along y-direction
     84    INTEGER(iwp) ::  k !< grid index along z-direction
    8385
    84        REAL(wp) ::  ukomp !< advection velocity along x-direction
    85        REAL(wp) ::  vkomp !< advection velocity along y-direction
    86        REAL(wp) ::  wkomp !< advection velocity along z-direction
     86    REAL(wp) ::  ukomp !< advection velocity along x-direction
     87    REAL(wp) ::  vkomp !< advection velocity along y-direction
     88    REAL(wp) ::  wkomp !< advection velocity along z-direction
    8789
    88        
    89        DO  i = nxlu, nxr
    90           DO  j = nys, nyn
    91              DO  k = nzb+1, nzt
     90
     91    DO  i = nxlu, nxr
     92       DO  j = nys, nyn
     93          DO  k = nzb+1, nzt
    9294!
    93 !--             x-direction
    94                 ukomp = u(k,j,i) - u_gtrans
    95                 IF ( ukomp > 0.0_wp )  THEN
    96                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    97                                          ( u(k,j,i) - u(k,j,i-1) ) * ddx
    98                 ELSE
    99                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    100                                           ( u(k,j,i+1) - u(k,j,i) ) * ddx
    101                 ENDIF
     95!--          x-direction
     96             ukomp = u(k,j,i) - u_gtrans
     97             IF ( ukomp > 0.0_wp )  THEN
     98                tend(k,j,i) = tend(k,j,i) - ukomp * ( u(k,j,i) - u(k,j,i-1) ) * ddx
     99             ELSE
     100                tend(k,j,i) = tend(k,j,i) - ukomp * ( u(k,j,i+1) - u(k,j,i) ) * ddx
     101             ENDIF
    102102!
    103 !--             y-direction
    104                 vkomp = 0.25_wp * ( v(k,j,i)   + v(k,j+1,i) +                  &
    105                                  v(k,j,i-1) + v(k,j+1,i-1) ) - v_gtrans
    106                 IF ( vkomp > 0.0_wp )  THEN
    107                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    108                                          ( u(k,j,i) - u(k,j-1,i) ) * ddy
    109                 ELSE
    110                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    111                                          ( u(k,j+1,i) - u(k,j,i) ) * ddy
    112                 ENDIF
     103!--          y-direction
     104             vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k,j,i-1) + v(k,j+1,i-1) ) - v_gtrans
     105             IF ( vkomp > 0.0_wp )  THEN
     106                tend(k,j,i) = tend(k,j,i) - vkomp * ( u(k,j,i) - u(k,j-1,i) ) * ddy
     107             ELSE
     108                tend(k,j,i) = tend(k,j,i) - vkomp * ( u(k,j+1,i) - u(k,j,i) ) * ddy
     109             ENDIF
    113110!
    114 !--             z-direction
    115                 wkomp = 0.25_wp * ( w(k,j,i)   + w(k-1,j,i) +                  &
    116                                  w(k,j,i-1) + w(k-1,j,i-1) )
    117                 IF ( wkomp > 0.0_wp )  THEN
    118                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    119                                          ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)
    120                 ELSE
    121                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    122                                          ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)
    123                 ENDIF
     111!--          z-direction
     112             wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j,i-1) + w(k-1,j,i-1) )
     113             IF ( wkomp > 0.0_wp )  THEN
     114                tend(k,j,i) = tend(k,j,i) - wkomp * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)
     115             ELSE
     116                tend(k,j,i) = tend(k,j,i) - wkomp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)
     117             ENDIF
    124118
    125              ENDDO
    126119          ENDDO
    127120       ENDDO
     121    ENDDO
    128122
    129     END SUBROUTINE advec_u_up
     123 END SUBROUTINE advec_u_up
    130124
    131125
    132 !------------------------------------------------------------------------------!
     126!--------------------------------------------------------------------------------------------------!
    133127! Description:
    134128! ------------
    135129!> Call for grid point i,j
    136 !------------------------------------------------------------------------------!
    137     SUBROUTINE advec_u_up_ij( i, j )
     130!--------------------------------------------------------------------------------------------------!
     131 SUBROUTINE advec_u_up_ij( i, j )
    138132
    139        USE arrays_3d,                                                          &
    140            ONLY:  ddzu, tend, u, v, w
     133    USE arrays_3d,                                                                                 &
     134        ONLY:  ddzu, tend, u, v, w
    141135
    142        USE control_parameters,                                                 &
    143            ONLY:  u_gtrans, v_gtrans
     136    USE control_parameters,                                                                        &
     137        ONLY:  u_gtrans, v_gtrans
    144138
    145        USE grid_variables,                                                     &
    146            ONLY:  ddx, ddy
     139    USE grid_variables,                                                                            &
     140        ONLY:  ddx, ddy
    147141
    148        USE indices,                                                            &
    149            ONLY:  nzb, nzt
     142    USE indices,                                                                                   &
     143        ONLY:  nzb, nzt
    150144
    151        USE kinds
     145    USE kinds
    152146
    153147
    154        IMPLICIT NONE
     148    IMPLICIT NONE
    155149
    156        INTEGER(iwp) ::  i !< grid index along x-direction
    157        INTEGER(iwp) ::  j !< grid index along y-direction
    158        INTEGER(iwp) ::  k !< grid index along z-direction
     150    INTEGER(iwp) ::  i !< grid index along x-direction
     151    INTEGER(iwp) ::  j !< grid index along y-direction
     152    INTEGER(iwp) ::  k !< grid index along z-direction
    159153
    160        REAL(wp) ::  ukomp !< advection velocity along x-direction
    161        REAL(wp) ::  vkomp !< advection velocity along y-direction
    162        REAL(wp) ::  wkomp !< advection velocity along z-direction
     154    REAL(wp) ::  ukomp !< advection velocity along x-direction
     155    REAL(wp) ::  vkomp !< advection velocity along y-direction
     156    REAL(wp) ::  wkomp !< advection velocity along z-direction
    163157
    164158
    165        DO  k = nzb+1, nzt
     159    DO  k = nzb+1, nzt
    166160!
    167 !--       x-direction
    168           ukomp = u(k,j,i) - u_gtrans
    169           IF ( ukomp > 0.0_wp )  THEN
    170              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    171                                          ( u(k,j,i) - u(k,j,i-1) ) * ddx
    172           ELSE
    173              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    174                                          ( u(k,j,i+1) - u(k,j,i) ) * ddx
    175           ENDIF
     161!--    x-direction
     162       ukomp = u(k,j,i) - u_gtrans
     163       IF ( ukomp > 0.0_wp )  THEN
     164          tend(k,j,i) = tend(k,j,i) - ukomp * ( u(k,j,i) - u(k,j,i-1) ) * ddx
     165       ELSE
     166          tend(k,j,i) = tend(k,j,i) - ukomp * ( u(k,j,i+1) - u(k,j,i) ) * ddx
     167       ENDIF
    176168!
    177 !--       y-direction
    178           vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k,j,i-1) + v(k,j+1,i-1) &
    179                          ) - v_gtrans
    180           IF ( vkomp > 0.0_wp )  THEN
    181              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    182                                          ( u(k,j,i) - u(k,j-1,i) ) * ddy
    183           ELSE
    184              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    185                                          ( u(k,j+1,i) - u(k,j,i) ) * ddy
    186           ENDIF
     169!--    y-direction
     170       vkomp = 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k,j,i-1) + v(k,j+1,i-1) ) - v_gtrans
     171       IF ( vkomp > 0.0_wp )  THEN
     172          tend(k,j,i) = tend(k,j,i) - vkomp * ( u(k,j,i) - u(k,j-1,i) ) * ddy
     173       ELSE
     174          tend(k,j,i) = tend(k,j,i) - vkomp * ( u(k,j+1,i) - u(k,j,i) ) * ddy
     175       ENDIF
    187176!
    188 !--       z-direction
    189           wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j,i-1) + w(k-1,j,i-1) )
    190           IF ( wkomp > 0.0_wp )  THEN
    191              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    192                                          ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)
    193           ELSE
    194              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    195                                          ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)
    196           ENDIF
     177!--    z-direction
     178       wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j,i-1) + w(k-1,j,i-1) )
     179       IF ( wkomp > 0.0_wp )  THEN
     180          tend(k,j,i) = tend(k,j,i) - wkomp * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k)
     181       ELSE
     182          tend(k,j,i) = tend(k,j,i) - wkomp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)
     183       ENDIF
    197184
    198        ENDDO
     185    ENDDO
    199186
    200     END SUBROUTINE advec_u_up_ij
     187 END SUBROUTINE advec_u_up_ij
    201188
    202189 END MODULE advec_u_up_mod
  • palm/trunk/SOURCE/advec_v_pw.f90

    r4360 r4488  
    11!> @file advec_v_pw.f90
    2 !------------------------------------------------------------------------------!
    3 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    43! This file is part of the PALM model system.
    54!
    6 ! PALM is free software: you can redistribute it and/or modify it under the
    7 ! terms of the GNU General Public License as published by the Free Software
    8 ! Foundation, either version 3 of the License, or (at your option) any later
    9 ! 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.
    108!
    11 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    12 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    13 ! 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.
    1412!
    15 ! You should have received a copy of the GNU General Public License along with
    16 ! 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/>.
    1715!
    1816! Copyright 1997-2020 Leibniz Universitaet Hannover
    19 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    2018!
    2119! Current revisions:
     
    2624! -----------------
    2725! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
    2829! Corrected "Former revisions" section
    2930!
     
    3839! ------------
    3940!> Advection term for v velocity-component using Piacsek and Williams.
    40 !> Vertical advection at the first grid point above the surface is done with
    41 !> normal centred differences, because otherwise no information from the surface
    42 !> would be communicated upwards due to w=0 at K=nzb.
    43 !------------------------------------------------------------------------------!
     41!> Vertical advection at the first grid point above the surface is done with normal centred
     42!> differences, because otherwise no information from the surface would be communicated upwards due
     43!> to w=0 at K=nzb.
     44!--------------------------------------------------------------------------------------------------!
    4445 MODULE advec_v_pw_mod
    4546 
     
    5657
    5758
    58 !------------------------------------------------------------------------------!
     59!--------------------------------------------------------------------------------------------------!
    5960! Description:
    6061! ------------
    6162!> Call for all grid points
    62 !------------------------------------------------------------------------------!
    63     SUBROUTINE advec_v_pw
     63!--------------------------------------------------------------------------------------------------!
     64 SUBROUTINE advec_v_pw
    6465
    65        USE arrays_3d,                                                          &
    66            ONLY:  ddzw, tend, u, v, w
     66    USE arrays_3d,                                                                                 &
     67        ONLY:  ddzw, tend, u, v, w
    6768
    68        USE control_parameters,                                                 &
    69            ONLY:  u_gtrans, v_gtrans
     69    USE control_parameters,                                                                        &
     70        ONLY:  u_gtrans, v_gtrans
    7071
    71        USE grid_variables,                                                     &
    72            ONLY:  ddx, ddy
     72    USE grid_variables,                                                                            &
     73        ONLY:  ddx, ddy
    7374
    74        USE indices,                                                            &
    75            ONLY:  nxl, nxr, nyn, nysv, nzb, nzt
     75    USE indices,                                                                                   &
     76        ONLY:  nxl, nxr, nyn, nysv, nzb, nzt
    7677
    77        USE kinds
     78    USE kinds
    7879
    7980
    80        IMPLICIT NONE
     81    IMPLICIT NONE
    8182
    82        INTEGER(iwp) ::  i !< grid index along x-direction
    83        INTEGER(iwp) ::  j !< grid index along y-direction
    84        INTEGER(iwp) ::  k !< grid index along z-direction
    85        
    86        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    87        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
    88  
     83    INTEGER(iwp) ::  i !< grid index along x-direction
     84    INTEGER(iwp) ::  j !< grid index along y-direction
     85    INTEGER(iwp) ::  k !< grid index along z-direction
     86   
     87    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     88    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
    8989
    90        gu = 2.0_wp * u_gtrans
    91        gv = 2.0_wp * v_gtrans
    92        DO  i = nxl, nxr
    93           DO  j = nysv, nyn
    94              DO  k = nzb+1, nzt
    95                 tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    96                          ( v(k,j,i+1) * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )     &
    97                          - v(k,j,i-1) * ( u(k,j-1,i) + u(k,j,i) - gu ) ) * ddx &
    98                        + ( v(k,j+1,i) * ( v(k,j+1,i) + v(k,j,i) - gv )         &
    99                          - v(k,j-1,i) * ( v(k,j,i) + v(k,j-1,i) - gv ) ) * ddy &
    100                        + ( v(k+1,j,i) * ( w(k,j-1,i) + w(k,j,i) )              &
    101                          - v(k-1,j,i) * ( w(k-1,j-1,i) + w(k-1,j,i) ) )        &
    102                                                                   * ddzw(k)    &
    103                                                       )
    104              ENDDO
     90
     91    gu = 2.0_wp * u_gtrans
     92    gv = 2.0_wp * v_gtrans
     93    DO  i = nxl, nxr
     94       DO  j = nysv, nyn
     95          DO  k = nzb+1, nzt
     96             tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                               &
     97                           ( v(k,j,i+1)   * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )                     &
     98                           - v(k,j,i-1)   * ( u(k,j-1,i)   + u(k,j,i)   - gu ) ) * ddx             &
     99                           + ( v(k,j+1,i) * ( v(k,j+1,i)   + v(k,j,i)   - gv )                     &
     100                           - v(k,j-1,i)   * ( v(k,j,i)     + v(k,j-1,i) - gv ) ) * ddy             &
     101                           + ( v(k+1,j,i) * ( w(k,j-1,i)   + w(k,j,i) )                            &
     102                           - v(k-1,j,i)   * ( w(k-1,j-1,i) + w(k-1,j,i) ) )      * ddzw(k)         &
     103                                                   )
    105104          ENDDO
    106105       ENDDO
     106    ENDDO
    107107
    108     END SUBROUTINE advec_v_pw
     108 END SUBROUTINE advec_v_pw
    109109
    110110
    111 !------------------------------------------------------------------------------!
     111!--------------------------------------------------------------------------------------------------!
    112112! Description:
    113113! ------------
    114114!> Call for grid point i,j
    115 !------------------------------------------------------------------------------!
    116     SUBROUTINE advec_v_pw_ij( i, j )
     115!--------------------------------------------------------------------------------------------------!
     116 SUBROUTINE advec_v_pw_ij( i, j )
    117117
    118        USE arrays_3d,                                                          &
    119            ONLY:  ddzw, tend, u, v, w
     118    USE arrays_3d,                                                                                 &
     119        ONLY:  ddzw, tend, u, v, w
    120120
    121        USE control_parameters,                                                 &
    122            ONLY:  u_gtrans, v_gtrans
     121    USE control_parameters,                                                                        &
     122        ONLY:  u_gtrans, v_gtrans
    123123
    124        USE grid_variables,                                                     &
    125            ONLY:  ddx, ddy
     124    USE grid_variables,                                                                            &
     125        ONLY:  ddx, ddy
    126126
    127        USE indices,                                                            &
    128            ONLY:  nzb, nzt
     127    USE indices,                                                                                   &
     128        ONLY:  nzb, nzt
    129129
    130        USE kinds
     130    USE kinds
    131131
    132132
    133        IMPLICIT NONE
     133    IMPLICIT NONE
    134134
    135        INTEGER(iwp) ::  i !< grid index along x-direction
    136        INTEGER(iwp) ::  j !< grid index along y-direction
    137        INTEGER(iwp) ::  k !< grid index along z-direction
    138        
    139        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    140        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
     135    INTEGER(iwp) ::  i !< grid index along x-direction
     136    INTEGER(iwp) ::  j !< grid index along y-direction
     137    INTEGER(iwp) ::  k !< grid index along z-direction
     138   
     139    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     140    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
    141141
    142142
    143        gu = 2.0_wp * u_gtrans
    144        gv = 2.0_wp * v_gtrans
    145        DO  k = nzb+1, nzt
    146           tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    147                          ( v(k,j,i+1) * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )     &
    148                          - v(k,j,i-1) * ( u(k,j-1,i) + u(k,j,i) - gu ) ) * ddx &
    149                        + ( v(k,j+1,i) * ( v(k,j+1,i) + v(k,j,i) - gv )         &
    150                          - v(k,j-1,i) * ( v(k,j,i) + v(k,j-1,i) - gv ) ) * ddy &
    151                        + ( v(k+1,j,i) * ( w(k,j-1,i) + w(k,j,i) )              &
    152                          - v(k-1,j,i) * ( w(k-1,j-1,i) + w(k-1,j,i) ) )        &
    153                                                                   * ddzw(k)    &
    154                                                 )
    155        ENDDO
    156  
    157     END SUBROUTINE advec_v_pw_ij
     143    gu = 2.0_wp * u_gtrans
     144    gv = 2.0_wp * v_gtrans
     145    DO  k = nzb+1, nzt
     146       tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                                     &
     147                      ( v(k,j,i+1)   * ( u(k,j-1,i+1) + u(k,j,i+1) - gu )                          &
     148                      - v(k,j,i-1)   * ( u(k,j-1,i)   + u(k,j,i)   - gu ) ) * ddx                  &
     149                      + ( v(k,j+1,i) * ( v(k,j+1,i)   + v(k,j,i)   - gv )                          &
     150                      - v(k,j-1,i)   * ( v(k,j,i)     + v(k,j-1,i) - gv ) ) * ddy                  &
     151                      + ( v(k+1,j,i) * ( w(k,j-1,i)   + w(k,j,i) )                                 &
     152                      - v(k-1,j,i)   * ( w(k-1,j-1,i) + w(k-1,j,i) ) )      * ddzw(k)              &
     153                                             )
     154    ENDDO
     155
     156 END SUBROUTINE advec_v_pw_ij
    158157
    159158 END MODULE advec_v_pw_mod
  • palm/trunk/SOURCE/advec_v_up.f90

    r4360 r4488  
    11!> @file advec_v_up.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
    2830!
     
    3941!> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0!
    4042!>       The same problem occurs for all topography boundaries!
    41 !------------------------------------------------------------------------------!
     43!--------------------------------------------------------------------------------------------------!
    4244 MODULE advec_v_up_mod
    4345 
     
    5456
    5557
    56 !------------------------------------------------------------------------------!
     58!--------------------------------------------------------------------------------------------------!
    5759! Description:
    5860! ------------
    5961!> Call for all grid points
    60 !------------------------------------------------------------------------------!
    61     SUBROUTINE advec_v_up
     62!--------------------------------------------------------------------------------------------------!
     63 SUBROUTINE advec_v_up
    6264
    63        USE arrays_3d,                                                          &
    64            ONLY:  ddzu, tend, u, v, w
     65    USE arrays_3d,                                                                                 &
     66        ONLY:  ddzu, tend, u, v, w
    6567
    66        USE control_parameters,                                                 &
    67            ONLY:  u_gtrans, v_gtrans
     68    USE control_parameters,                                                                        &
     69        ONLY:  u_gtrans, v_gtrans
    6870
    69        USE grid_variables,                                                     &
    70            ONLY:  ddx, ddy
     71    USE grid_variables,                                                                            &
     72        ONLY:  ddx, ddy
    7173
    72        USE indices,                                                            &
    73            ONLY:  nxl, nxr, nyn, nysv, nzb, nzt
     74    USE indices,                                                                                   &
     75        ONLY:  nxl, nxr, nyn, nysv, nzb, nzt
    7476
    75        USE kinds
     77    USE kinds
    7678
    7779
    78        IMPLICIT NONE
     80    IMPLICIT NONE
    7981
    80        INTEGER(iwp) ::  i !< grid index along x-direction
    81        INTEGER(iwp) ::  j !< grid index along y-direction
    82        INTEGER(iwp) ::  k !< grid index along z-direction
     82    INTEGER(iwp) ::  i !< grid index along x-direction
     83    INTEGER(iwp) ::  j !< grid index along y-direction
     84    INTEGER(iwp) ::  k !< grid index along z-direction
    8385
    84        REAL(wp) ::  ukomp !< advection velocity along x-direction
    85        REAL(wp) ::  vkomp !< advection velocity along y-direction
    86        REAL(wp) ::  wkomp !< advection velocity along z-direction       
     86    REAL(wp) ::  ukomp !< advection velocity along x-direction
     87    REAL(wp) ::  vkomp !< advection velocity along y-direction
     88    REAL(wp) ::  wkomp !< advection velocity along z-direction       
    8789
    8890
    89        DO  i = nxl, nxr
    90           DO  j = nysv, nyn
    91              DO  k = nzb+1, nzt
     91    DO  i = nxl, nxr
     92       DO  j = nysv, nyn
     93          DO  k = nzb+1, nzt
    9294!
    93 !--             x-direction
    94                 ukomp = 0.25_wp * ( u(k,j,i)   + u(k,j-1,i) +                  &
    95                                  u(k,j,i+1) + u(k,j-1,i+1) ) - u_gtrans
    96                 IF ( ukomp > 0.0_wp )  THEN
    97                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    98                                          ( v(k,j,i) - v(k,j,i-1) ) * ddx
    99                 ELSE
    100                    tend(k,j,i) = tend(k,j,i) - ukomp *                         &
    101                                          ( v(k,j,i+1) - v(k,j,i) ) * ddx
    102                 ENDIF
     95!--          x-direction
     96             ukomp = 0.25_wp * ( u(k,j,i) + u(k,j-1,i) + u(k,j,i+1) + u(k,j-1,i+1) ) - u_gtrans
     97             IF ( ukomp > 0.0_wp )  THEN
     98                tend(k,j,i) = tend(k,j,i) - ukomp * ( v(k,j,i) - v(k,j,i-1) ) * ddx
     99             ELSE
     100                tend(k,j,i) = tend(k,j,i) - ukomp * ( v(k,j,i+1) - v(k,j,i) ) * ddx
     101             ENDIF
    103102!
    104 !--             y-direction
    105                 vkomp = v(k,j,i) - v_gtrans
    106                 IF ( vkomp > 0.0_wp )  THEN
    107                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    108                                          ( v(k,j,i) - v(k,j-1,i) ) * ddy
    109                 ELSE
    110                    tend(k,j,i) = tend(k,j,i) - vkomp *                         &
    111                                          ( v(k,j+1,i) - v(k,j,i) ) * ddy
    112                 ENDIF
     103!--          y-direction
     104             vkomp = v(k,j,i) - v_gtrans
     105             IF ( vkomp > 0.0_wp )  THEN
     106                tend(k,j,i) = tend(k,j,i) - vkomp * ( v(k,j,i) - v(k,j-1,i) ) * ddy
     107             ELSE
     108                tend(k,j,i) = tend(k,j,i) - vkomp * ( v(k,j+1,i) - v(k,j,i) ) * ddy
     109             ENDIF
    113110!
    114 !--             z-direction
    115                 wkomp = 0.25_wp * ( w(k,j,i)  + w(k-1,j,i) +                   &
    116                                  w(k,j-1,i) + w(k-1,j-1,i) )
    117                 IF ( wkomp > 0.0_wp )  THEN
    118                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    119                                          ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)
    120                 ELSE
    121                    tend(k,j,i) = tend(k,j,i) - wkomp *                         &
    122                                          ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)
    123                 ENDIF
     111!--          z-direction
     112             wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )
     113             IF ( wkomp > 0.0_wp )  THEN
     114                tend(k,j,i) = tend(k,j,i) - wkomp * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)
     115             ELSE
     116                tend(k,j,i) = tend(k,j,i) - wkomp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)
     117             ENDIF
    124118
    125              ENDDO
    126119          ENDDO
    127120       ENDDO
     121    ENDDO
    128122
    129     END SUBROUTINE advec_v_up
     123 END SUBROUTINE advec_v_up
    130124
    131125
    132 !------------------------------------------------------------------------------!
     126!--------------------------------------------------------------------------------------------------!
    133127! Description:
    134128! ------------
    135129!> Call for grid point i,j
    136 !------------------------------------------------------------------------------!
    137     SUBROUTINE advec_v_up_ij( i, j )
     130!--------------------------------------------------------------------------------------------------!
     131 SUBROUTINE advec_v_up_ij( i, j )
    138132
    139        USE arrays_3d,                                                          &
    140            ONLY:  ddzu, tend, u, v, w
     133    USE arrays_3d,                                                                                 &
     134        ONLY:  ddzu, tend, u, v, w
    141135
    142        USE control_parameters,                                                 &
    143            ONLY:  u_gtrans, v_gtrans
     136    USE control_parameters,                                                                        &
     137        ONLY:  u_gtrans, v_gtrans
    144138
    145        USE grid_variables,                                                     &
    146            ONLY:  ddx, ddy
     139    USE grid_variables,                                                                            &
     140        ONLY:  ddx, ddy
    147141
    148        USE indices,                                                            &
    149            ONLY:  nzb, nzt
     142    USE indices,                                                                                   &
     143        ONLY:  nzb, nzt
    150144
    151        USE kinds
     145    USE kinds
    152146
    153147
    154        IMPLICIT NONE
     148    IMPLICIT NONE
    155149
    156        INTEGER(iwp) ::  i !< grid index along x-direction
    157        INTEGER(iwp) ::  j !< grid index along y-direction
    158        INTEGER(iwp) ::  k !< grid index along z-direction
     150    INTEGER(iwp) ::  i !< grid index along x-direction
     151    INTEGER(iwp) ::  j !< grid index along y-direction
     152    INTEGER(iwp) ::  k !< grid index along z-direction
    159153
    160        REAL(wp) ::  ukomp !< advection velocity along x-direction
    161        REAL(wp) ::  vkomp !< advection velocity along y-direction
    162        REAL(wp) ::  wkomp !< advection velocity along z-direction
     154    REAL(wp) ::  ukomp !< advection velocity along x-direction
     155    REAL(wp) ::  vkomp !< advection velocity along y-direction
     156    REAL(wp) ::  wkomp !< advection velocity along z-direction
    163157
    164158
    165        DO  k = nzb+1, nzt
     159    DO  k = nzb+1, nzt
    166160!
    167 !--       x-direction
    168           ukomp = 0.25_wp * ( u(k,j,i) + u(k,j-1,i) + u(k,j,i+1) + u(k,j-1,i+1) &
    169                          ) - u_gtrans
    170           IF ( ukomp > 0.0_wp )  THEN
    171              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    172                                          ( v(k,j,i) - v(k,j,i-1) ) * ddx
    173           ELSE
    174              tend(k,j,i) = tend(k,j,i) - ukomp *                               &
    175                                          ( v(k,j,i+1) - v(k,j,i) ) * ddx
    176           ENDIF
     161!--    x-direction
     162       ukomp = 0.25_wp * ( u(k,j,i) + u(k,j-1,i) + u(k,j,i+1) + u(k,j-1,i+1) ) - u_gtrans
     163       IF ( ukomp > 0.0_wp )  THEN
     164          tend(k,j,i) = tend(k,j,i) - ukomp * ( v(k,j,i) - v(k,j,i-1) ) * ddx
     165       ELSE
     166          tend(k,j,i) = tend(k,j,i) - ukomp * ( v(k,j,i+1) - v(k,j,i) ) * ddx
     167       ENDIF
    177168!
    178 !--       y-direction
    179           vkomp = v(k,j,i) - v_gtrans
    180           IF ( vkomp > 0.0_wp )  THEN
    181              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    182                                          ( v(k,j,i) - v(k,j-1,i) ) * ddy
    183           ELSE
    184              tend(k,j,i) = tend(k,j,i) - vkomp *                               &
    185                                          ( v(k,j+1,i) - v(k,j,i) ) * ddy
    186           ENDIF
     169!--    y-direction
     170       vkomp = v(k,j,i) - v_gtrans
     171       IF ( vkomp > 0.0_wp )  THEN
     172          tend(k,j,i) = tend(k,j,i) - vkomp * ( v(k,j,i) - v(k,j-1,i) ) * ddy
     173       ELSE
     174          tend(k,j,i) = tend(k,j,i) - vkomp * ( v(k,j+1,i) - v(k,j,i) ) * ddy
     175       ENDIF
    187176!
    188 !--       z-direction
    189           wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )
    190           IF ( wkomp > 0.0_wp )  THEN
    191              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    192                                          ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)
    193           ELSE
    194              tend(k,j,i) = tend(k,j,i) - wkomp *                               &
    195                                          ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)
    196           ENDIF
     177!--    z-direction
     178       wkomp = 0.25_wp * ( w(k,j,i) + w(k-1,j,i) + w(k,j-1,i) + w(k-1,j-1,i) )
     179       IF ( wkomp > 0.0_wp )  THEN
     180          tend(k,j,i) = tend(k,j,i) - wkomp * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)
     181       ELSE
     182          tend(k,j,i) = tend(k,j,i) - wkomp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)
     183       ENDIF
    197184
    198        ENDDO
     185    ENDDO
    199186
    200     END SUBROUTINE advec_v_up_ij
     187 END SUBROUTINE advec_v_up_ij
    201188
    202189 END MODULE advec_v_up_mod
  • palm/trunk/SOURCE/advec_w_pw.f90

    r4360 r4488  
    11!> @file advec_w_pw.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
    2830!
     
    3739! ------------
    3840!> Advection term for w velocity-component using Piacsek and Williams.
    39 !> Vertical advection at the first grid point above the surface is done with
    40 !> normal centred differences, because otherwise no information from the surface
    41 !> would be communicated upwards due to w=0 at k=nzb.
    42 !------------------------------------------------------------------------------!
     41!> Vertical advection at the first grid point above the surface is done with normal centred
     42!> differences, because otherwise no information from the surface would be communicated upwards due
     43!> to w=0 at k=nzb.
     44!--------------------------------------------------------------------------------------------------!
    4345 MODULE advec_w_pw_mod
    4446 
     
    5557
    5658
    57 !------------------------------------------------------------------------------!
     59!--------------------------------------------------------------------------------------------------!
    5860! Description:
    5961! ------------
    6062!> Call for all grid points
    61 !------------------------------------------------------------------------------!
    62     SUBROUTINE advec_w_pw
     63!--------------------------------------------------------------------------------------------------!
     64 SUBROUTINE advec_w_pw
    6365
    64        USE arrays_3d,                                                          &
    65            ONLY:  ddzu, tend, u, v, w
     66    USE arrays_3d,                                                                                 &
     67        ONLY:  ddzu, tend, u, v, w
    6668
    67        USE control_parameters,                                                 &
    68            ONLY:  u_gtrans, v_gtrans
     69    USE control_parameters,                                                                        &
     70        ONLY:  u_gtrans, v_gtrans
    6971
    70        USE grid_variables,                                                     &
    71            ONLY:  ddx, ddy
     72    USE grid_variables,                                                                            &
     73        ONLY:  ddx, ddy
    7274
    73        USE indices,                                                            &
    74            ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     75    USE indices,                                                                                   &
     76        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
    7577
    76        USE kinds
     78    USE kinds
    7779
    7880
    79        IMPLICIT NONE
     81    IMPLICIT NONE
    8082
    81        INTEGER(iwp) ::  i !< grid index along x-direction
    82        INTEGER(iwp) ::  j !< grid index along y-direction
    83        INTEGER(iwp) ::  k !< grid index along z-direction
    84        
    85        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    86        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
     83    INTEGER(iwp) ::  i !< grid index along x-direction
     84    INTEGER(iwp) ::  j !< grid index along y-direction
     85    INTEGER(iwp) ::  k !< grid index along z-direction
     86   
     87    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     88    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
    8789
    88  
    89        gu = 2.0_wp * u_gtrans
    90        gv = 2.0_wp * v_gtrans
    91        DO  i = nxl, nxr
    92           DO  j = nys, nyn
    93              DO  k = nzb+1, nzt
    94                 tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                        &
    95                          ( w(k,j,i+1) * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )     &
    96                          - w(k,j,i-1) * ( u(k+1,j,i) + u(k,j,i) - gu ) ) * ddx &
    97                        + ( w(k,j+1,i) * ( v(k+1,j+1,i) + v(k,j+1,i) - gv )     &
    98                          - w(k,j-1,i) * ( v(k+1,j,i) + v(k,j,i) - gv ) ) * ddy &
    99                        + ( w(k+1,j,i) * ( w(k+1,j,i) + w(k,j,i) )              &
    100                          - w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i) ) )            &
    101                                                                   * ddzu(k+1)  &
    102                                                       )
    103              ENDDO
     90
     91    gu = 2.0_wp * u_gtrans
     92    gv = 2.0_wp * v_gtrans
     93    DO  i = nxl, nxr
     94       DO  j = nys, nyn
     95          DO  k = nzb+1, nzt
     96             tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                               &
     97                           ( w(k,j,i+1)   * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )                     &
     98                           - w(k,j,i-1)   * ( u(k+1,j,i)   + u(k,j,i) - gu ) ) * ddx               &
     99                           + ( w(k,j+1,i) * ( v(k+1,j+1,i) + v(k,j+1,i) - gv )                     &
     100                           - w(k,j-1,i)   * ( v(k+1,j,i)   + v(k,j,i) - gv ) ) * ddy               &
     101                           + ( w(k+1,j,i) * ( w(k+1,j,i)   + w(k,j,i) )                            &
     102                           - w(k-1,j,i)   * ( w(k,j,i)     + w(k-1,j,i) ) )    * ddzu(k+1)         &
     103                                                   )
    104104          ENDDO
    105105       ENDDO
     106    ENDDO
    106107
    107     END SUBROUTINE advec_w_pw
     108 END SUBROUTINE advec_w_pw
    108109
    109110
    110 !------------------------------------------------------------------------------!
     111!--------------------------------------------------------------------------------------------------!
    111112! Description:
    112113! ------------
    113114!> Call for grid point i,j
    114 !------------------------------------------------------------------------------!
    115     SUBROUTINE advec_w_pw_ij( i, j )
     115!--------------------------------------------------------------------------------------------------!
     116 SUBROUTINE advec_w_pw_ij( i, j )
    116117
    117        USE arrays_3d,                                                          &
    118            ONLY:  ddzu, tend, u, v, w
     118    USE arrays_3d,                                                                                 &
     119        ONLY:  ddzu, tend, u, v, w
    119120
    120        USE control_parameters,                                                 &
    121            ONLY:  u_gtrans, v_gtrans
     121    USE control_parameters,                                                                        &
     122        ONLY:  u_gtrans, v_gtrans
    122123
    123        USE grid_variables,                                                     &
    124            ONLY:  ddx, ddy
     124    USE grid_variables,                                                                            &
     125        ONLY:  ddx, ddy
    125126
    126        USE indices,                                                            &
    127            ONLY:  nzb, nzt
     127    USE indices,                                                                                   &
     128        ONLY:  nzb, nzt
    128129
    129        USE kinds
     130    USE kinds
    130131
    131132
    132        IMPLICIT NONE
     133    IMPLICIT NONE
    133134
    134        INTEGER(iwp) ::  i !< grid index along x-direction
    135        INTEGER(iwp) ::  j !< grid index along y-direction
    136        INTEGER(iwp) ::  k !< grid index along z-direction
    137        
    138        REAL(wp)    ::  gu !< Galilei-transformation velocity along x
    139        REAL(wp)    ::  gv !< Galilei-transformation velocity along y
     135    INTEGER(iwp) ::  i !< grid index along x-direction
     136    INTEGER(iwp) ::  j !< grid index along y-direction
     137    INTEGER(iwp) ::  k !< grid index along z-direction
     138   
     139    REAL(wp)     ::  gu !< Galilei-transformation velocity along x
     140    REAL(wp)     ::  gv !< Galilei-transformation velocity along y
    140141
    141        gu = 2.0_wp * u_gtrans
    142        gv = 2.0_wp * v_gtrans
    143        DO  k = nzb+1, nzt
    144           tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                              &
    145                          ( w(k,j,i+1) * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )     &
    146                          - w(k,j,i-1) * ( u(k+1,j,i) + u(k,j,i) - gu ) ) * ddx &
    147                        + ( w(k,j+1,i) * ( v(k+1,j+1,i) + v(k,j+1,i) - gv )     &
    148                          - w(k,j-1,i) * ( v(k+1,j,i) + v(k,j,i) - gv ) ) * ddy &
    149                        + ( w(k+1,j,i) * ( w(k+1,j,i) + w(k,j,i) )              &
    150                          - w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i) ) )            &
    151                                                                   * ddzu(k+1)  &
    152                                                 )
    153        ENDDO
    154     END SUBROUTINE advec_w_pw_ij
     142    gu = 2.0_wp * u_gtrans
     143    gv = 2.0_wp * v_gtrans
     144    DO  k = nzb+1, nzt
     145       tend(k,j,i) = tend(k,j,i) - 0.25_wp * (                                                     &
     146                      ( w(k,j,i+1)   * ( u(k+1,j,i+1) + u(k,j,i+1) - gu )                          &
     147                      - w(k,j,i-1)   * ( u(k+1,j,i)   + u(k,j,i) - gu ) ) * ddx                    &
     148                      + ( w(k,j+1,i) * ( v(k+1,j+1,i) + v(k,j+1,i) - gv )                          &
     149                      - w(k,j-1,i)   * ( v(k+1,j,i)   + v(k,j,i) - gv ) ) * ddy                    &
     150                      + ( w(k+1,j,i) * ( w(k+1,j,i)   + w(k,j,i) )                                 &
     151                      - w(k-1,j,i)   * ( w(k,j,i)     + w(k-1,j,i) ) )    * ddzu(k+1)              &
     152                                             )
     153    ENDDO
     154 END SUBROUTINE advec_w_pw_ij
    155155
    156156 END MODULE advec_w_pw_mod
Note: See TracChangeset for help on using the changeset viewer.