source: palm/trunk/SOURCE/advec_s_up.f90 @ 4281

Last change on this file since 4281 was 4182, checked in by scharf, 22 months ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 7.4 KB
RevLine 
[1873]1!> @file advec_s_up.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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/>.
16!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1354]22!
[2233]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: advec_s_up.f90 4182 2019-08-22 15:20:23Z schwenkel $
[4182]27! Corrected "Former revisions" section
28!
29! 3927 2019-04-23 13:24:29Z raasch
[3927]30! pointer attribute removed from scalar 3d-array for performance reasons
31!
32! 3665 2019-01-10 08:28:24Z raasch
[3665]33! unused variables removed
34!
35! 3655 2019-01-07 16:51:22Z knoop
[3636]36! nopointer option removed
[1321]37!
[4182]38! Revision 1.1  1997/08/29 08:54:33  raasch
39! Initial revision
40!
41!
[1]42! Description:
43! ------------
[1682]44!> Advection term for scalar quantities using the Upstream scheme.
45!> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0!
46!>       The same problem occurs for all topography boundaries!
[1]47!------------------------------------------------------------------------------!
[1682]48 MODULE advec_s_up_mod
49 
[1]50
51    PRIVATE
52    PUBLIC advec_s_up
53
54    INTERFACE advec_s_up
55       MODULE PROCEDURE advec_s_up
56       MODULE PROCEDURE advec_s_up_ij
57    END INTERFACE advec_s_up
58
59 CONTAINS
60
61
62!------------------------------------------------------------------------------!
[1682]63! Description:
64! ------------
65!> Call for all grid points
[1]66!------------------------------------------------------------------------------!
67    SUBROUTINE advec_s_up( sk )
68
[1320]69       USE arrays_3d,                                                          &
70           ONLY:  ddzu, tend, u, v, w
[1]71
[1320]72       USE control_parameters,                                                 &
73           ONLY:  u_gtrans, v_gtrans
74
75       USE grid_variables,                                                     &
76           ONLY:  ddx, ddy
77
78       USE indices,                                                            &
[3927]79           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
[1320]80
81       USE kinds
82
83
[1]84       IMPLICIT NONE
85
[3547]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
[1]89
[3547]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
[3636]93
[3927]94       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
[1]95
96
97       DO  i = nxl, nxr
98          DO  j = nys, nyn
[2232]99             DO  k = nzb+1, nzt
[1]100!
101!--             x-direction
[1353]102                ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
103                IF ( ukomp > 0.0_wp )  THEN
[1320]104                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
[3538]105                                         ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
[1]106                ELSE
[1320]107                   tend(k,j,i) = tend(k,j,i) - ukomp *                         &
[3538]108                                         ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
[1]109                ENDIF
110!
111!--             y-direction
[1353]112                vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
113                IF ( vkomp > 0.0_wp )  THEN
[1320]114                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
[3538]115                                         ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
[1]116                ELSE
[1320]117                   tend(k,j,i) = tend(k,j,i) - vkomp *                         &
[3538]118                                         ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
[1]119                ENDIF
120!
121!--             z-direction
[1353]122                wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
123                IF ( wkomp > 0.0_wp )  THEN
[1320]124                   tend(k,j,i) = tend(k,j,i) - wkomp *                         &
[3538]125                                         ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
[1]126                ELSE
[1320]127                   tend(k,j,i) = tend(k,j,i) - wkomp *                         &
[3538]128                                         ( sk(k+1,j,i)-sk(k,j,i) ) * ddzu(k+1)
[1]129                ENDIF
130
131             ENDDO
132          ENDDO
133       ENDDO
134
135    END SUBROUTINE advec_s_up
136
137
138!------------------------------------------------------------------------------!
[1682]139! Description:
140! ------------
141!> Call for grid point i,j
[1]142!------------------------------------------------------------------------------!
143    SUBROUTINE advec_s_up_ij( i, j, sk )
144
[1320]145       USE arrays_3d,                                                          &
146           ONLY:  ddzu, tend, u, v, w
[1]147
[1320]148       USE control_parameters,                                                 &
149           ONLY:  u_gtrans, v_gtrans
150
151       USE grid_variables,                                                     &
152           ONLY:  ddx, ddy
153
154       USE indices,                                                            &
[3927]155           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
[1320]156
157       USE kinds
158
159
[1]160       IMPLICIT NONE
161
[3547]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
[1]165
[3547]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
[1320]169       
[3927]170       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< treated scalar
[1]171
172
[2232]173       DO  k = nzb+1, nzt
[1]174!
175!--       x-direction
[1353]176          ukomp = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - u_gtrans
177          IF ( ukomp > 0.0_wp )  THEN
[1320]178             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
[3538]179                                         ( sk(k,j,i) - sk(k,j,i-1) ) * ddx
[1]180          ELSE
[1320]181             tend(k,j,i) = tend(k,j,i) - ukomp *                               &
[3538]182                                         ( sk(k,j,i+1) - sk(k,j,i) ) * ddx
[1]183          ENDIF
184!
185!--       y-direction
[1353]186          vkomp = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - v_gtrans
187          IF ( vkomp > 0.0_wp )  THEN
[1320]188             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
[3538]189                                         ( sk(k,j,i) - sk(k,j-1,i) ) * ddy
[1]190          ELSE
[1320]191             tend(k,j,i) = tend(k,j,i) - vkomp *                               &
[3538]192                                         ( sk(k,j+1,i) - sk(k,j,i) ) * ddy
[1]193          ENDIF
194!
195!--       z-direction
[1353]196          wkomp = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
197          IF ( wkomp > 0.0_wp )  THEN
[1320]198             tend(k,j,i) = tend(k,j,i) - wkomp *                               &
[3538]199                                         ( sk(k,j,i) - sk(k-1,j,i) ) * ddzu(k)
[1]200          ELSE
[1320]201             tend(k,j,i) = tend(k,j,i) - wkomp *                               &
[3538]202                                         ( sk(k+1,j,i)-sk(k,j,i) ) * ddzu(k+1)
[1]203          ENDIF
204
205       ENDDO
206
207    END SUBROUTINE advec_s_up_ij
208
209 END MODULE advec_s_up_mod
Note: See TracBrowser for help on using the repository browser.