Changeset 4488 for palm/trunk
- Timestamp:
- Apr 3, 2020 11:34:29 AM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 8 edited
-
advec_s_bc.f90 (modified) (70 diffs)
-
advec_s_pw.f90 (modified) (4 diffs)
-
advec_s_up.f90 (modified) (4 diffs)
-
advec_u_pw.f90 (modified) (4 diffs)
-
advec_u_up.f90 (modified) (4 diffs)
-
advec_v_pw.f90 (modified) (4 diffs)
-
advec_v_up.f90 (modified) (4 diffs)
-
advec_w_pw.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_bc.f90
r4429 r4488 1 1 !> @file advec_s_bc.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 ! ------------------------------------------------------------------------------!17 ! -------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4429 2020-02-27 15:24:30Z raasch 27 29 ! bugfix: cpp-directives added for serial mode 28 ! 30 ! 29 31 ! 4360 2020-01-07 11:25:50Z suehring 30 32 ! Corrected "Former revisions" section 31 ! 33 ! 32 34 ! 3761 2019-02-25 15:31:42Z raasch 33 35 ! unused variables removed 34 ! 36 ! 35 37 ! 3655 2019-01-07 16:51:22Z knoop 36 38 ! nopointer option removed … … 45 47 !> Computation in individual steps for each of the three dimensions. 46 48 !> 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. 52 53 !> NOTE: This is a provisional, non-optimised version! 53 !------------------------------------------------------------------------------ !54 MODULE advec_s_bc_mod55 54 !--------------------------------------------------------------------------------------------------! 55 MODULE advec_s_bc_mod 56 56 57 57 58 PRIVATE … … 64 65 CONTAINS 65 66 66 !------------------------------------------------------------------------------ !67 !--------------------------------------------------------------------------------------------------! 67 68 ! Description: 68 69 ! ------------ 69 70 !> @todo Missing subroutine description. 70 !------------------------------------------------------------------------------ !71 !--------------------------------------------------------------------------------------------------! 71 72 SUBROUTINE advec_s_bc( sk, sk_char ) 72 73 73 USE advection, &74 USE advection, & 74 75 ONLY: aex, bex, dex, eex 75 76 76 USE arrays_3d, &77 USE arrays_3d, & 77 78 ONLY: d, ddzw, dzu, dzw, tend, u, v, w 78 79 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, & 85 85 ONLY: cpu_log, log_point_s 86 86 87 USE grid_variables, &87 USE grid_variables, & 88 88 ONLY: ddx, ddy 89 89 90 USE indices, &90 USE indices, & 91 91 ONLY: nx, nxl, nxr, nyn, nys, nzb, nzt 92 92 … … 95 95 USE pegrid 96 96 97 USE statistics, &97 USE statistics, & 98 98 ONLY: rmask, statistic_regions, sums_wsts_bc_l 99 99 … … 112 112 INTEGER(iwp) :: type_xz_2 !< 113 113 #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 !< 142 141 143 142 REAL(wp) :: fmax(2) !< 144 143 REAL(wp) :: fmax_l(2) !< 145 146 REAL(wp), DIMENSION(:,:,:), POINTER :: sk147 144 148 145 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a0 !< … … 161 158 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1 !< 162 159 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sw !< 163 160 164 161 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sk_p !< 162 163 REAL(wp), DIMENSION(:,:,:), POINTER :: sk !< 164 165 165 166 166 ! … … 197 197 ! 198 198 !-- 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, & 201 201 comm2d, status, ierr ) 202 202 ! 203 203 !-- 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, & 206 206 comm2d, status, ierr ) 207 207 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) … … 217 217 218 218 ! 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 total221 !-- d omain must be adjusted by the temperature difference between this distance219 !-- 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 222 222 IF ( sloping_surface .AND. sk_char == 'pt' ) THEN 223 223 IF ( nxl == 0 ) THEN … … 241 241 DO j = nys, nyn 242 242 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) ) 244 244 denomi = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 245 245 fmax_l(1) = MAX( fmax_l(1) , numera ) … … 259 259 ! 260 260 !-- 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) & 269 268 ) 270 269 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 271 270 272 271 ! 273 !-- Initialise point of time measuring of the exponential portion (this would 274 !-- not work if donelocally within the loop)272 !-- Initialise point of time measuring of the exponential portion (this would not work if done 273 !-- locally within the loop) 275 274 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' ) 276 275 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) 277 276 278 !279 !-- Outer loop of all j280 277 DO j = nys, nyn 281 278 … … 285 282 DO k = nzb+1, nzt 286 283 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)) * f1920292 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) & 294 291 ) * 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 298 296 ENDDO 299 297 ENDDO … … 301 299 ! 302 300 !-- Fluxes using the Bott scheme 303 !-- *VOCL LOOP,UNROLL(2)304 301 DO i = nxl, nxr 305 302 DO k = nzb+1, nzt … … 308 305 cipf = 1.0_wp - 2.0_wp * cip 309 306 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 ) & 312 309 + 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 ) & 315 312 + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 316 313 ip = MAX( ip, 0.0_wp ) … … 323 320 cipf = 1.0_wp - 2.0_wp * cip 324 321 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 ) & 327 324 + 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 ) & 330 327 + a2(k,i) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 331 328 ip = MAX( ip, 0.0_wp ) … … 358 355 DO i = nxl-1, nxr+1 359 356 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 ) 362 358 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0_wp 363 359 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 ) 366 361 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0_wp 367 362 … … 369 364 t2 = 0.35_wp 370 365 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 ) & 377 369 ) sw(k,i) = 1.0_wp 378 370 ENDDO … … 385 377 DO k = nzb+1, nzt 386 378 387 !-- *VOCL STMT,IF(10)388 379 IF ( sw(k,i) == 1.0_wp ) THEN 389 380 snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1) … … 397 388 cip = MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 398 389 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 ) & 404 395 ) 405 396 IF ( sterm == 0.0001_wp ) ippe(k,i) = sk_p(k,j,i) * cip … … 416 407 cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 417 408 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 ) & 423 414 ) 424 415 IF ( sterm == 0.0001_wp ) imme(k,i) = sk_p(k,j,i) * cim … … 426 417 ENDIF 427 418 428 !-- *VOCL STMT,IF(10)429 419 IF ( sw(k,i+1) == 1.0_wp ) THEN 430 420 snenn = sk_p(k,j,i) - sk_p(k,j,i+2) … … 438 428 cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) 439 429 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 ) & 445 435 ) 446 436 IF ( sterm == 0.0001_wp ) impe(k,i) = sk_p(k,j,i+1) * cim … … 448 438 ENDIF 449 439 450 !-- *VOCL STMT,IF(10)451 440 IF ( sw(k,i-1) == 1.0_wp ) THEN 452 441 snenn = sk_p(k,j,i) - sk_p(k,j,i-2) … … 460 449 cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) 461 450 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 ) & 467 456 ) 468 457 IF ( sterm == 0.0001_wp ) ipme(k,i) = sk_p(k,j,i-1) * cip … … 478 467 DO i = nxl, nxr 479 468 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) 484 473 tendcy = fplus - fminus 485 474 ! … … 490 479 !-- Density correction because of possible remaining divergences 491 480 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 ) 494 482 d(k,j,i) = d_new 495 483 ENDDO … … 500 488 ! 501 489 !-- 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 ) 504 491 505 492 … … 508 495 #if defined( __parallel ) 509 496 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, & 511 498 type_xz_2, ierr ) 512 499 CALL MPI_TYPE_COMMIT( type_xz_2, ierr ) … … 514 501 !-- Send front boundary, receive rear boundary 515 502 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, & 518 505 comm2d, status, ierr ) 519 506 ! 520 507 !-- 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, & 523 510 comm2d, status, ierr ) 524 511 CALL MPI_TYPE_FREE( type_xz_2, ierr ) … … 561 548 ! 562 549 !-- 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) & 571 557 ) 572 558 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 573 559 574 !575 !-- Outer loop of all i576 560 DO i = nxl, nxr 577 561 … … 581 565 DO k = nzb+1, nzt 582 566 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)) * f1920588 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) & 590 574 ) * 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 594 579 ENDDO 595 580 ENDDO … … 597 582 ! 598 583 !-- Fluxes using the Bott scheme 599 !-- *VOCL LOOP,UNROLL(2)600 584 DO j = nys, nyn 601 585 DO k = nzb+1, nzt … … 604 588 cipf = 1.0_wp - 2.0_wp * cip 605 589 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 ) & 608 592 + 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 ) & 611 595 + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 612 596 ip = MAX( ip, 0.0_wp ) … … 619 603 cipf = 1.0_wp - 2.0_wp * cip 620 604 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 ) & 623 607 + 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 ) & 626 610 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 627 611 ip = MAX( ip, 0.0_wp ) … … 654 638 DO j = nys-1, nyn+1 655 639 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) ) / & 657 641 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 658 642 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 659 643 660 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / &644 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 661 645 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 662 646 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp … … 666 650 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 667 651 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 ) & 673 655 ) sw(k,j) = 1.0_wp 674 656 ENDDO … … 681 663 DO k = nzb+1, nzt 682 664 683 !-- *VOCL STMT,IF(10)684 665 IF ( sw(k,j) == 1.0_wp ) THEN 685 666 snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i) … … 693 674 cip = MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 694 675 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 ) & 700 681 ) 701 682 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip … … 712 693 cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 713 694 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 ) & 719 700 ) 720 701 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim … … 722 703 ENDIF 723 704 724 !-- *VOCL STMT,IF(10)725 705 IF ( sw(k,j+1) == 1.0_wp ) THEN 726 706 snenn = sk_p(k,j,i) - sk_p(k,j+2,i) … … 734 714 cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) 735 715 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 ) & 741 721 ) 742 722 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k,j+1,i) * cim … … 744 724 ENDIF 745 725 746 !-- *VOCL STMT,IF(10)747 726 IF ( sw(k,j-1) == 1.0_wp ) THEN 748 727 snenn = sk_p(k,j,i) - sk_p(k,j-2,i) … … 756 735 cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) 757 736 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 ) & 763 742 ) 764 743 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k,j-1,i) * cip … … 774 753 DO j = nys, nyn 775 754 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) 780 759 tendcy = fplus - fminus 781 760 ! … … 786 765 !-- Density correction because of possible remaining divergences 787 766 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 791 769 ENDDO 792 770 ENDDO 793 771 794 772 ENDDO ! End of the advection in y-direction 773 774 795 775 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 796 776 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' ) … … 798 778 ! 799 779 !-- 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 ) 802 781 803 782 … … 959 938 960 939 ! 961 !-- Cloud water content boundary condition at the bottom boundary: 962 !-- Dirichlet (fixed surfacerain water content).940 !-- Cloud water content boundary condition at the bottom boundary: Dirichlet (fixed surface 941 !-- rain water content). 963 942 DO i = nxl, nxr 964 943 DO j = nys, nyn … … 981 960 982 961 ! 983 !-- Rain water content boundary condition at the bottom boundary: 984 !-- Dirichlet (fixed surfacerain water content).962 !-- Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface 963 !-- rain water content). 985 964 DO i = nxl, nxr 986 965 DO j = nys, nyn … … 1003 982 1004 983 ! 1005 !-- Cloud drop concentration boundary condition at the bottom boundary: 1006 !-- Dirichlet (fixedsurface cloud drop concentration).984 !-- Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed 985 !-- surface cloud drop concentration). 1007 986 DO i = nxl, nxr 1008 987 DO j = nys, nyn … … 1025 1004 1026 1005 ! 1027 !-- Rain drop concentration boundary condition at the bottom boundary: 1028 !-- Dirichlet (fixedsurface rain drop concentration).1006 !-- Rain drop concentration boundary condition at the bottom boundary: Dirichlet (fixed 1007 !-- surface rain drop concentration). 1029 1008 DO i = nxl, nxr 1030 1009 DO j = nys, nyn … … 1060 1039 ELSE 1061 1040 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, '"' 1064 1042 CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 ) 1065 1043 1066 1044 ENDIF 1067 1045 … … 1090 1068 ! 1091 1069 !-- 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) & 1100 1075 ) 1101 1076 imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp 1102 1077 1103 !1104 !-- Outer loop of all i1105 1078 DO i = nxl, nxr 1106 1079 … … 1110 1083 DO k = nzb, nzt+1 1111 1084 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)) * f19201117 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) & 1119 1092 ) * 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 1123 1097 ENDDO 1124 1098 ENDDO … … 1126 1100 ! 1127 1101 !-- Fluxes using the Bott scheme 1128 !-- *VOCL LOOP,UNROLL(2)1129 1102 DO j = nys, nyn 1130 1103 DO k = nzb+1, nzt … … 1133 1106 cipf = 1.0_wp - 2.0_wp * cip 1134 1107 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 ) & 1137 1110 + 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 ) & 1140 1113 + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1141 1114 ip = MAX( ip, 0.0_wp ) … … 1148 1121 cipf = 1.0_wp - 2.0_wp * cip 1149 1122 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 ) & 1152 1125 + 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 ) & 1155 1128 + a2(k,j) * f24 * ( 1.0_wp - cimf*cimf*cimf ) 1156 1129 ip = MAX( ip, 0.0_wp ) … … 1183 1156 DO j = nys, nyn 1184 1157 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) ) / & 1186 1159 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp ) 1187 1160 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0_wp 1188 1161 1189 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / &1162 m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) / & 1190 1163 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp ) 1191 1164 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0_wp … … 1195 1168 IF ( m1(k,j) == -1.0_wp ) t2 = 0.12_wp 1196 1169 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 ) & 1202 1173 ) sw(k,j) = 1.0_wp 1203 1174 ENDDO … … 1210 1181 DO k = nzb+1, nzt 1211 1182 1212 !-- *VOCL STMT,IF(10)1213 1183 IF ( sw(k,j) == 1.0_wp ) THEN 1214 1184 snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i) … … 1222 1192 cip = MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1223 1193 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 ) & 1229 1199 ) 1230 1200 IF ( sterm == 0.0001_wp ) ippe(k,j) = sk_p(k,j,i) * cip … … 1241 1211 cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1242 1212 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 ) & 1248 1218 ) 1249 1219 IF ( sterm == 0.0001_wp ) imme(k,j) = sk_p(k,j,i) * cim … … 1251 1221 ENDIF 1252 1222 1253 !-- *VOCL STMT,IF(10)1254 1223 IF ( sw(k+1,j) == 1.0_wp ) THEN 1255 1224 snenn = sk_p(k,j,i) - sk_p(k+2,j,i) … … 1263 1232 cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) ) 1264 1233 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 ) & 1270 1239 ) 1271 1240 IF ( sterm == 0.0001_wp ) impe(k,j) = sk_p(k+1,j,i) * cim … … 1273 1242 ENDIF 1274 1243 1275 !-- *VOCL STMT,IF(10)1276 1244 IF ( sw(k-1,j) == 1.0_wp ) THEN 1277 1245 snenn = sk_p(k,j,i) - sk_p(k-2,j,i) … … 1285 1253 cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) ) 1286 1254 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 ) & 1292 1260 ) 1293 1261 IF ( sterm == 0.0001_wp ) ipme(k,j) = sk_p(k-1,j,i) * cip … … 1303 1271 DO j = nys, nyn 1304 1272 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) 1309 1277 tendcy = fplus - fminus 1310 1278 ! … … 1315 1283 !-- Density correction because of possible remaining divergences 1316 1284 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 ) 1319 1286 ! 1320 1287 !-- Store heat flux for subsequent statistics output. 1321 !-- array m1 is here used as temporary storage1288 !-- Array m1 is here used as temporary storage 1322 1289 m1(k,j) = fplus / dt_3d * dzw(k) 1323 1290 ENDDO … … 1325 1292 1326 1293 ! 1327 !-- Sum up heat flux in order to o rder to obtain horizontal averages1294 !-- Sum up heat flux in order to obtain horizontal averages 1328 1295 IF ( sk_char == 'pt' ) THEN 1329 1296 DO sr = 0, statistic_regions 1330 1297 DO j = nys, nyn 1331 1298 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) 1334 1300 ENDDO 1335 1301 ENDDO … … 1337 1303 ENDIF 1338 1304 1339 ENDDO ! End of the advection in z-direction 1305 ENDDO 1306 ! 1307 !-- End of the advection in z-direction 1308 1340 1309 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) 1341 1310 CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' ) … … 1343 1312 ! 1344 1313 !-- 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 ) 1347 1315 1348 1316 ! -
palm/trunk/SOURCE/advec_s_pw.f90
r4360 r4488 1 1 !> @file advec_s_pw.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 ! 30 ! 29 31 ! 3927 2019-04-23 13:24:29Z raasch 30 32 ! pointer attribute removed from scalar 3d-array for performance reasons 31 ! 33 ! 32 34 ! 3665 2019-01-10 08:28:24Z raasch 33 35 ! unused variables removed 34 ! 36 ! 35 37 ! 3655 2019-01-07 16:51:22Z knoop 36 38 ! nopointer option removed … … 42 44 ! Description: 43 45 ! ------------ 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)). 49 50 !> NOTE: at the first grid point above the surface computation still takes place! 50 !------------------------------------------------------------------------------ !51 !--------------------------------------------------------------------------------------------------! 51 52 MODULE advec_s_pw_mod 52 53 53 54 54 55 PRIVATE … … 59 60 MODULE PROCEDURE advec_s_pw_ij 60 61 END INTERFACE 61 62 62 63 CONTAINS 63 64 64 65 65 !------------------------------------------------------------------------------ !66 !--------------------------------------------------------------------------------------------------! 66 67 ! Description: 67 68 ! ------------ 68 69 !> Call for all grid points 69 !------------------------------------------------------------------------------ !70 70 !--------------------------------------------------------------------------------------------------! 71 SUBROUTINE advec_s_pw( sk ) 71 72 72 USE arrays_3d,&73 73 USE arrays_3d, & 74 ONLY: dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w 74 75 75 USE control_parameters,&76 76 USE control_parameters, & 77 ONLY: u_gtrans, v_gtrans 77 78 78 USE grid_variables,&79 79 USE grid_variables, & 80 ONLY: ddx, ddy 80 81 81 USE indices,&82 82 USE indices, & 83 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 83 84 84 85 USE kinds 85 86 86 87 87 88 IMPLICIT NONE 88 89 89 90 91 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 92 93 93 94 94 REAL(wp) :: gu !< local additional advective velocity 95 REAL(wp) :: gv !< local additional advective velocity 95 96 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 98 98 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 102 103 103 104 ! 104 !-- 105 106 105 !-- Galilean transformation + Stokes drift velocity (ocean only) 106 gu = u_gtrans - u_stokes_zu(k) 107 gv = v_gtrans - v_stokes_zu(k) 107 108 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 ) 126 126 ENDDO 127 127 ENDDO 128 ENDDO 128 129 129 130 END SUBROUTINE advec_s_pw 130 131 131 132 132 !------------------------------------------------------------------------------ !133 !--------------------------------------------------------------------------------------------------! 133 134 ! Description: 134 135 ! ------------ 135 136 !> Call for grid point i,j 136 !------------------------------------------------------------------------------ !137 137 !--------------------------------------------------------------------------------------------------! 138 SUBROUTINE advec_s_pw_ij( i, j, sk ) 138 139 139 USE arrays_3d,&140 140 USE arrays_3d, & 141 ONLY: dd2zu, tend, u, u_stokes_zu, v, v_stokes_zu, w 141 142 142 USE control_parameters,&143 143 USE control_parameters, & 144 ONLY: u_gtrans, v_gtrans 144 145 145 USE grid_variables,&146 146 USE grid_variables, & 147 ONLY: ddx, ddy 147 148 148 USE indices,&149 149 USE indices, & 150 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt 150 151 151 152 USE kinds 152 153 153 154 154 155 IMPLICIT NONE 155 156 156 157 158 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 159 160 160 161 161 REAL(wp) :: gu !< local additional advective velocity 162 REAL(wp) :: gv !< local additional advective velocity 162 163 163 164 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk 164 165 165 166 166 167 DO k = nzb+1, nzt 167 168 168 169 ! 169 !-- 170 171 170 !-- Galilean transformation + Stokes drift velocity (ocean only) 171 gu = u_gtrans - u_stokes_zu(k) 172 gv = v_gtrans - v_stokes_zu(k) 172 173 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 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 191 192 192 193 END SUBROUTINE advec_s_pw_ij 193 194 194 195 END MODULE advec_s_pw_mod -
palm/trunk/SOURCE/advec_s_up.f90
r4360 r4488 1 1 !> @file advec_s_up.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 ! 30 ! 29 31 ! 3927 2019-04-23 13:24:29Z raasch 30 32 ! pointer attribute removed from scalar 3d-array for performance reasons 31 ! 33 ! 32 34 ! 3665 2019-01-10 08:28:24Z raasch 33 35 ! unused variables removed 34 ! 36 ! 35 37 ! 3655 2019-01-07 16:51:22Z knoop 36 38 ! nopointer option removed … … 45 47 !> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0! 46 48 !> The same problem occurs for all topography boundaries! 47 !------------------------------------------------------------------------------ !49 !--------------------------------------------------------------------------------------------------! 48 50 MODULE advec_s_up_mod 49 51 50 52 51 53 PRIVATE … … 60 62 61 63 62 !------------------------------------------------------------------------------ !64 !--------------------------------------------------------------------------------------------------! 63 65 ! Description: 64 66 ! ------------ 65 67 !> Call for all grid points 66 !------------------------------------------------------------------------------ !67 68 !--------------------------------------------------------------------------------------------------! 69 SUBROUTINE advec_s_up( sk ) 68 70 69 USE arrays_3d,&70 71 USE arrays_3d, & 72 ONLY: ddzu, tend, u, v, w 71 73 72 USE control_parameters,&73 74 USE control_parameters, & 75 ONLY: u_gtrans, v_gtrans 74 76 75 USE grid_variables,&76 77 USE grid_variables, & 78 ONLY: ddx, ddy 77 79 78 USE indices,&79 80 USE indices, & 81 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 80 82 81 83 USE kinds 82 84 83 85 84 86 IMPLICIT NONE 85 87 86 87 88 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 89 91 90 REAL(wp) :: ukomp !< advection velocity along x-direction91 REAL(wp) :: vkomp !< advection velocity along y-direction92 REAL(wp) :: wkomp !< advection velocity along z-direction92 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 93 95 94 96 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< treated scalar 95 97 96 98 97 98 99 99 DO i = nxl, nxr 100 DO j = nys, nyn 101 DO k = nzb+1, nzt 100 102 ! 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 110 110 ! 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 120 118 ! 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 130 126 131 ENDDO132 127 ENDDO 133 128 ENDDO 129 ENDDO 134 130 135 131 END SUBROUTINE advec_s_up 136 132 137 133 138 !------------------------------------------------------------------------------ !134 !--------------------------------------------------------------------------------------------------! 139 135 ! Description: 140 136 ! ------------ 141 137 !> Call for grid point i,j 142 !------------------------------------------------------------------------------ !143 138 !--------------------------------------------------------------------------------------------------! 139 SUBROUTINE advec_s_up_ij( i, j, sk ) 144 140 145 USE arrays_3d,&146 141 USE arrays_3d, & 142 ONLY: ddzu, tend, u, v, w 147 143 148 USE control_parameters,&149 144 USE control_parameters, & 145 ONLY: u_gtrans, v_gtrans 150 146 151 USE grid_variables,&152 147 USE grid_variables, & 148 ONLY: ddx, ddy 153 149 154 USE indices,&155 150 USE indices, & 151 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt 156 152 157 153 USE kinds 158 154 159 155 160 156 IMPLICIT NONE 161 157 162 163 164 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 165 161 166 REAL(wp) :: ukomp !< advection velocity along x-direction167 REAL(wp) :: vkomp !< advection velocity along y-direction168 REAL(wp) :: wkomp !< advection velocity along z-direction169 170 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 171 167 172 168 173 169 DO k = nzb+1, nzt 174 170 ! 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 184 178 ! 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 194 186 ! 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 204 194 205 195 ENDDO 206 196 207 197 END SUBROUTINE advec_s_up_ij 208 198 209 199 END MODULE advec_s_up_mod -
palm/trunk/SOURCE/advec_u_pw.f90
r4360 r4488 1 1 !> @file advec_u_pw.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 ! 30 ! 29 31 ! 3655 2019-01-07 16:51:22Z knoop 30 32 ! variables documented … … 37 39 ! ------------ 38 40 !> 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 surface41 !> would be communicated upwards dueto 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 !--------------------------------------------------------------------------------------------------! 43 45 MODULE advec_u_pw_mod 44 46 45 47 46 48 PRIVATE … … 51 53 MODULE PROCEDURE advec_u_pw_ij 52 54 END INTERFACE advec_u_pw 53 55 54 56 CONTAINS 55 57 56 58 57 !------------------------------------------------------------------------------ !59 !--------------------------------------------------------------------------------------------------! 58 60 ! Description: 59 61 ! ------------ 60 62 !> Call for all grid points 61 !------------------------------------------------------------------------------ !62 63 !--------------------------------------------------------------------------------------------------! 64 SUBROUTINE advec_u_pw 63 65 64 USE arrays_3d,&65 66 USE arrays_3d, & 67 ONLY: ddzw, tend, u, v, w 66 68 67 USE control_parameters,&68 69 USE control_parameters, & 70 ONLY: u_gtrans, v_gtrans 69 71 70 USE grid_variables,&71 72 USE grid_variables, & 73 ONLY: ddx, ddy 72 74 73 USE indices,&74 75 USE indices, & 76 ONLY: nxlu, nxr, nyn, nys, nzb, nzt 75 77 76 78 USE kinds 77 79 78 80 79 81 IMPLICIT NONE 80 82 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 ) 103 103 ENDDO 104 104 ENDDO 105 ENDDO 105 106 106 107 END SUBROUTINE advec_u_pw 107 108 108 109 109 !------------------------------------------------------------------------------ !110 !--------------------------------------------------------------------------------------------------! 110 111 ! Description: 111 112 ! ------------ 112 113 !> Call for grid point i,j 113 !------------------------------------------------------------------------------ !114 114 !--------------------------------------------------------------------------------------------------! 115 SUBROUTINE advec_u_pw_ij( i, j ) 115 116 116 USE arrays_3d,&117 117 USE arrays_3d, & 118 ONLY: ddzw, tend, u, v, w 118 119 119 USE control_parameters,&120 120 USE control_parameters, & 121 ONLY: u_gtrans, v_gtrans 121 122 122 USE grid_variables,&123 123 USE grid_variables, & 124 ONLY: ddx, ddy 124 125 125 USE indices,&126 126 USE indices, & 127 ONLY: nzb, nzt 127 128 128 129 USE kinds 129 130 130 131 131 132 IMPLICIT NONE 132 133 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 139 137 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 153 140 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 155 155 156 156 END MODULE advec_u_pw_mod -
palm/trunk/SOURCE/advec_u_up.f90
r4360 r4488 1 1 !> @file advec_u_up.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 ! 30 ! 29 31 ! 3655 2019-01-07 16:51:22Z knoop 30 32 ! variables documented … … 39 41 !> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0! 40 42 !> The same problem occurs for all topography boundaries! 41 !------------------------------------------------------------------------------ !43 !--------------------------------------------------------------------------------------------------! 42 44 MODULE advec_u_up_mod 43 45 44 46 45 47 PRIVATE … … 54 56 55 57 56 !------------------------------------------------------------------------------ !58 !--------------------------------------------------------------------------------------------------! 57 59 ! Description: 58 60 ! ------------ 59 61 !> Call for all grid points 60 !------------------------------------------------------------------------------ !61 62 !--------------------------------------------------------------------------------------------------! 63 SUBROUTINE advec_u_up 62 64 63 USE arrays_3d,&64 65 USE arrays_3d, & 66 ONLY: ddzu, tend, u, v, w 65 67 66 USE control_parameters,&67 68 USE control_parameters, & 69 ONLY: u_gtrans, v_gtrans 68 70 69 USE grid_variables,&70 71 USE grid_variables, & 72 ONLY: ddx, ddy 71 73 72 USE indices,&73 74 USE indices, & 75 ONLY: nxlu, nxr, nyn, nys, nzb, nzt 74 76 75 77 USE kinds 76 78 77 79 78 80 IMPLICIT NONE 79 81 80 81 82 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 83 85 84 REAL(wp) :: ukomp !< advection velocity along x-direction85 REAL(wp) :: vkomp !< advection velocity along y-direction86 REAL(wp) :: wkomp !< advection velocity along z-direction86 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 87 89 88 89 90 91 90 91 DO i = nxlu, nxr 92 DO j = nys, nyn 93 DO k = nzb+1, nzt 92 94 ! 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 102 102 ! 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 113 110 ! 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 124 118 125 ENDDO126 119 ENDDO 127 120 ENDDO 121 ENDDO 128 122 129 123 END SUBROUTINE advec_u_up 130 124 131 125 132 !------------------------------------------------------------------------------ !126 !--------------------------------------------------------------------------------------------------! 133 127 ! Description: 134 128 ! ------------ 135 129 !> Call for grid point i,j 136 !------------------------------------------------------------------------------ !137 130 !--------------------------------------------------------------------------------------------------! 131 SUBROUTINE advec_u_up_ij( i, j ) 138 132 139 USE arrays_3d,&140 133 USE arrays_3d, & 134 ONLY: ddzu, tend, u, v, w 141 135 142 USE control_parameters,&143 136 USE control_parameters, & 137 ONLY: u_gtrans, v_gtrans 144 138 145 USE grid_variables,&146 139 USE grid_variables, & 140 ONLY: ddx, ddy 147 141 148 USE indices,&149 142 USE indices, & 143 ONLY: nzb, nzt 150 144 151 145 USE kinds 152 146 153 147 154 148 IMPLICIT NONE 155 149 156 157 158 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 159 153 160 REAL(wp) :: ukomp !< advection velocity along x-direction161 REAL(wp) :: vkomp !< advection velocity along y-direction162 REAL(wp) :: wkomp !< advection velocity along z-direction154 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 163 157 164 158 165 159 DO k = nzb+1, nzt 166 160 ! 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 176 168 ! 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 187 176 ! 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 197 184 198 185 ENDDO 199 186 200 187 END SUBROUTINE advec_u_up_ij 201 188 202 189 END MODULE advec_u_up_mod -
palm/trunk/SOURCE/advec_v_pw.f90
r4360 r4488 1 1 !> @file advec_v_pw.f90 2 !------------------------------------------------------------------------------! 3 !------------------------------------------------------------------------------! 2 !--------------------------------------------------------------------------------------------------! 4 3 ! This file is part of the PALM model system. 5 4 ! 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. 10 8 ! 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 FOR13 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 14 12 ! 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/>. 17 15 ! 18 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 19 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 20 18 ! 21 19 ! Current revisions: … … 26 24 ! ----------------- 27 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 28 29 ! Corrected "Former revisions" section 29 30 ! … … 38 39 ! ------------ 39 40 !> 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 surface42 !> would be communicated upwards dueto 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 !--------------------------------------------------------------------------------------------------! 44 45 MODULE advec_v_pw_mod 45 46 … … 56 57 57 58 58 !------------------------------------------------------------------------------ !59 !--------------------------------------------------------------------------------------------------! 59 60 ! Description: 60 61 ! ------------ 61 62 !> Call for all grid points 62 !------------------------------------------------------------------------------ !63 63 !--------------------------------------------------------------------------------------------------! 64 SUBROUTINE advec_v_pw 64 65 65 USE arrays_3d,&66 66 USE arrays_3d, & 67 ONLY: ddzw, tend, u, v, w 67 68 68 USE control_parameters,&69 69 USE control_parameters, & 70 ONLY: u_gtrans, v_gtrans 70 71 71 USE grid_variables,&72 72 USE grid_variables, & 73 ONLY: ddx, ddy 73 74 74 USE indices,&75 75 USE indices, & 76 ONLY: nxl, nxr, nyn, nysv, nzb, nzt 76 77 77 78 USE kinds 78 79 79 80 80 81 IMPLICIT NONE 81 82 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 89 89 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 ) 105 104 ENDDO 106 105 ENDDO 106 ENDDO 107 107 108 108 END SUBROUTINE advec_v_pw 109 109 110 110 111 !------------------------------------------------------------------------------ !111 !--------------------------------------------------------------------------------------------------! 112 112 ! Description: 113 113 ! ------------ 114 114 !> Call for grid point i,j 115 !------------------------------------------------------------------------------ !116 115 !--------------------------------------------------------------------------------------------------! 116 SUBROUTINE advec_v_pw_ij( i, j ) 117 117 118 USE arrays_3d,&119 118 USE arrays_3d, & 119 ONLY: ddzw, tend, u, v, w 120 120 121 USE control_parameters,&122 121 USE control_parameters, & 122 ONLY: u_gtrans, v_gtrans 123 123 124 USE grid_variables,&125 124 USE grid_variables, & 125 ONLY: ddx, ddy 126 126 127 USE indices,&128 127 USE indices, & 128 ONLY: nzb, nzt 129 129 130 130 USE kinds 131 131 132 132 133 133 IMPLICIT NONE 134 134 135 136 137 138 139 REAL(wp):: gu !< Galilei-transformation velocity along x140 REAL(wp):: gv !< Galilei-transformation velocity along y135 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 141 141 142 142 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 158 157 159 158 END MODULE advec_v_pw_mod -
palm/trunk/SOURCE/advec_v_up.f90
r4360 r4488 1 1 !> @file advec_v_up.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 30 ! … … 39 41 !> NOTE: vertical advection at k=1 still has wrong grid spacing for w>0! 40 42 !> The same problem occurs for all topography boundaries! 41 !------------------------------------------------------------------------------ !43 !--------------------------------------------------------------------------------------------------! 42 44 MODULE advec_v_up_mod 43 45 … … 54 56 55 57 56 !------------------------------------------------------------------------------ !58 !--------------------------------------------------------------------------------------------------! 57 59 ! Description: 58 60 ! ------------ 59 61 !> Call for all grid points 60 !------------------------------------------------------------------------------ !61 62 !--------------------------------------------------------------------------------------------------! 63 SUBROUTINE advec_v_up 62 64 63 USE arrays_3d,&64 65 USE arrays_3d, & 66 ONLY: ddzu, tend, u, v, w 65 67 66 USE control_parameters,&67 68 USE control_parameters, & 69 ONLY: u_gtrans, v_gtrans 68 70 69 USE grid_variables,&70 71 USE grid_variables, & 72 ONLY: ddx, ddy 71 73 72 USE indices,&73 74 USE indices, & 75 ONLY: nxl, nxr, nyn, nysv, nzb, nzt 74 76 75 77 USE kinds 76 78 77 79 78 80 IMPLICIT NONE 79 81 80 81 82 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 83 85 84 85 86 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 87 89 88 90 89 90 91 91 DO i = nxl, nxr 92 DO j = nysv, nyn 93 DO k = nzb+1, nzt 92 94 ! 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 103 102 ! 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 113 110 ! 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 124 118 125 ENDDO126 119 ENDDO 127 120 ENDDO 121 ENDDO 128 122 129 123 END SUBROUTINE advec_v_up 130 124 131 125 132 !------------------------------------------------------------------------------ !126 !--------------------------------------------------------------------------------------------------! 133 127 ! Description: 134 128 ! ------------ 135 129 !> Call for grid point i,j 136 !------------------------------------------------------------------------------ !137 130 !--------------------------------------------------------------------------------------------------! 131 SUBROUTINE advec_v_up_ij( i, j ) 138 132 139 USE arrays_3d,&140 133 USE arrays_3d, & 134 ONLY: ddzu, tend, u, v, w 141 135 142 USE control_parameters,&143 136 USE control_parameters, & 137 ONLY: u_gtrans, v_gtrans 144 138 145 USE grid_variables,&146 139 USE grid_variables, & 140 ONLY: ddx, ddy 147 141 148 USE indices,&149 142 USE indices, & 143 ONLY: nzb, nzt 150 144 151 145 USE kinds 152 146 153 147 154 148 IMPLICIT NONE 155 149 156 157 158 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 159 153 160 161 162 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 163 157 164 158 165 159 DO k = nzb+1, nzt 166 160 ! 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 177 168 ! 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 187 176 ! 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 197 184 198 185 ENDDO 199 186 200 187 END SUBROUTINE advec_v_up_ij 201 188 202 189 END MODULE advec_v_up_mod -
palm/trunk/SOURCE/advec_w_pw.f90
r4360 r4488 1 1 !> @file advec_w_pw.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 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. 9 8 ! 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 FOR12 ! A PARTICULAR PURPOSE. See the GNU GeneralPublic 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. 13 12 ! 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/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4360 2020-01-07 11:25:50Z suehring 27 29 ! Corrected "Former revisions" section 28 30 ! … … 37 39 ! ------------ 38 40 !> 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 surface41 !> would be communicated upwards dueto 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 !--------------------------------------------------------------------------------------------------! 43 45 MODULE advec_w_pw_mod 44 46 … … 55 57 56 58 57 !------------------------------------------------------------------------------ !59 !--------------------------------------------------------------------------------------------------! 58 60 ! Description: 59 61 ! ------------ 60 62 !> Call for all grid points 61 !------------------------------------------------------------------------------ !62 63 !--------------------------------------------------------------------------------------------------! 64 SUBROUTINE advec_w_pw 63 65 64 USE arrays_3d,&65 66 USE arrays_3d, & 67 ONLY: ddzu, tend, u, v, w 66 68 67 USE control_parameters,&68 69 USE control_parameters, & 70 ONLY: u_gtrans, v_gtrans 69 71 70 USE grid_variables,&71 72 USE grid_variables, & 73 ONLY: ddx, ddy 72 74 73 USE indices,&74 75 USE indices, & 76 ONLY: nxl, nxr, nyn, nys, nzb, nzt 75 77 76 78 USE kinds 77 79 78 80 79 81 IMPLICIT NONE 80 82 81 82 83 84 85 REAL(wp):: gu !< Galilei-transformation velocity along x86 REAL(wp):: gv !< Galilei-transformation velocity along y83 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 87 89 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 ) 104 104 ENDDO 105 105 ENDDO 106 ENDDO 106 107 107 108 END SUBROUTINE advec_w_pw 108 109 109 110 110 !------------------------------------------------------------------------------ !111 !--------------------------------------------------------------------------------------------------! 111 112 ! Description: 112 113 ! ------------ 113 114 !> Call for grid point i,j 114 !------------------------------------------------------------------------------ !115 115 !--------------------------------------------------------------------------------------------------! 116 SUBROUTINE advec_w_pw_ij( i, j ) 116 117 117 USE arrays_3d,&118 118 USE arrays_3d, & 119 ONLY: ddzu, tend, u, v, w 119 120 120 USE control_parameters,&121 121 USE control_parameters, & 122 ONLY: u_gtrans, v_gtrans 122 123 123 USE grid_variables,&124 124 USE grid_variables, & 125 ONLY: ddx, ddy 125 126 126 USE indices,&127 127 USE indices, & 128 ONLY: nzb, nzt 128 129 129 130 USE kinds 130 131 131 132 132 133 IMPLICIT NONE 133 134 134 135 136 137 138 REAL(wp):: gu !< Galilei-transformation velocity along x139 REAL(wp):: gv !< Galilei-transformation velocity along y135 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 140 141 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 155 155 156 156 END MODULE advec_w_pw_mod
Note: See TracChangeset
for help on using the changeset viewer.
