Changeset 1320 for palm/trunk/SOURCE/advec_s_bc.f90
- Timestamp:
- Mar 20, 2014 8:40:49 AM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_s_bc.f90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_s_bc.f90
r1319 r1320 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! ONLY-attribute added to USE-statements, 23 ! kind-parameters added to all INTEGER and REAL declaration statements, 24 ! kinds are defined in new module kinds, 25 ! revision history before 2012 removed, 26 ! comment fields (!:) to be used for variable explanations added to 27 ! all variable declaration statements 23 28 ! 24 29 ! Former revisions: … … 37 42 ! 1010 2012-09-20 07:59:54Z raasch 38 43 ! cpp switch __nopointer added for pointer free version 39 !40 ! 622 2010-12-10 08:08:13Z raasch41 ! optional barriers included in order to speed up collective operations42 !43 ! 247 2009-02-27 14:01:30Z heinze44 ! Output of messages replaced by message handling routine45 !46 ! 216 2008-11-25 07:12:43Z raasch47 ! Neumann boundary condition at k=nzb is explicitly set for better reading,48 ! although this has been already done in boundary_conds49 !50 ! 97 2007-06-21 08:23:15Z raasch51 ! Advection of salinity included52 ! Bugfix: Error in boundary condition for TKE removed53 !54 ! 63 2007-03-13 03:52:49Z raasch55 ! Calculation extended for gridpoint nzt56 !57 ! RCS Log replace by Id keyword, revision history cleaned up58 !59 ! Revision 1.22 2006/02/23 09:42:08 raasch60 ! anz renamed ngp61 44 ! 62 45 ! Revision 1.1 1997/08/29 08:53:46 raasch … … 77 60 !------------------------------------------------------------------------------! 78 61 79 USE advection 80 USE arrays_3d 81 USE control_parameters 82 USE cpulog 83 USE grid_variables 84 USE indices 62 USE advection, & 63 ONLY: aex, bex, dex, eex 64 65 USE arrays_3d, & 66 ONLY: d, ddzw, dzu, dzw, tend, u, v, w 67 68 USE control_parameters, & 69 ONLY: dt_3d, bc_pt_t_val, bc_q_t_val, ibc_pt_b, ibc_pt_t, ibc_q_t, & 70 message_string, pt_slope_offset, sloping_surface, u_gtrans, & 71 v_gtrans 72 73 USE cpulog, & 74 ONLY: cpu_log, log_point_s 75 76 USE grid_variables, & 77 ONLY: ddx, ddy 78 79 USE indices, & 80 ONLY: nx, nxl, nxr, nyn, nys, nzb, nzt 81 82 USE kinds 83 85 84 USE pegrid 86 USE statistics 85 86 USE statistics, & 87 ONLY: rmask, statistic_regions, sums_wsts_bc_l 88 87 89 88 90 IMPLICIT NONE 89 91 90 CHARACTER (LEN=*) :: sk_char 91 92 INTEGER :: i, ix, j, k, ngp, sr, type_xz_2 93 94 REAL :: cim, cimf, cip, cipf, d_new, fminus, fplus, f2, f4, f8, & 95 f12, f24, f48, f1920, im, ip, m2, m3, nenner, snenn, sterm, & 96 tendenz, t1, t2, zaehler 97 REAL :: fmax(2), fmax_l(2) 92 CHARACTER (LEN=*) :: sk_char !: 93 94 INTEGER(iwp) :: i !: 95 INTEGER(iwp) :: ix !: 96 INTEGER(iwp) :: j !: 97 INTEGER(iwp) :: k !: 98 INTEGER(iwp) :: ngp !: 99 INTEGER(iwp) :: sr !: 100 INTEGER(iwp) :: type_xz_2 !: 101 102 REAL(wp) :: cim !: 103 REAL(wp) :: cimf !: 104 REAL(wp) :: cip !: 105 REAL(wp) :: cipf !: 106 REAL(wp) :: d_new !: 107 REAL(wp) :: denomi !: denominator 108 REAL(wp) :: fminus !: 109 REAL(wp) :: fplus !: 110 REAL(wp) :: f2 !: 111 REAL(wp) :: f4 !: 112 REAL(wp) :: f8 !: 113 REAL(wp) :: f12 !: 114 REAL(wp) :: f24 !: 115 REAL(wp) :: f48 !: 116 REAL(wp) :: f1920 !: 117 REAL(wp) :: im !: 118 REAL(wp) :: ip !: 119 REAL(wp) :: m2 !: 120 REAL(wp) :: m3 !: 121 REAL(wp) :: numera !: numerator 122 REAL(wp) :: snenn !: 123 REAL(wp) :: sterm !: 124 REAL(wp) :: tendcy !: 125 REAL(wp) :: t1 !: 126 REAL(wp) :: t2 !: 127 128 REAL(wp) :: fmax(2) !: 129 REAL(wp) :: fmax_l(2) !: 130 98 131 #if defined( __nopointer ) 99 REAL , DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk132 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: 100 133 #else 101 REAL , DIMENSION(:,:,:), POINTER :: sk134 REAL(wp), DIMENSION(:,:,:), POINTER :: sk 102 135 #endif 103 136 104 REAL, DIMENSION(:,:), ALLOCATABLE :: a0, a1, a12, a2, a22, immb, imme, & 105 impb, impe, ipmb, ipme, ippb, ippe 106 REAL, DIMENSION(:,:,:), ALLOCATABLE :: sk_p 137 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a0 !: 138 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a1 !: 139 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a12 !: 140 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a2 !: 141 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: a22 !: 142 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: immb !: 143 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: imme !: 144 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impb !: 145 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: impe !: 146 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipmb !: 147 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ipme !: 148 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippb !: 149 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ippe !: 150 151 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: sk_p !: 107 152 108 153 #if defined( __nec ) 109 REAL (kind=4) :: m1n, m1z !Wichtig: Division110 REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw154 REAL(sp) :: m1n, m1z !Wichtig: Division !: 155 REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !: 111 156 #else 112 REAL :: m1n, m1z113 REAL , DIMENSION(:,:), ALLOCATABLE :: m1, sw157 REAL(wp) :: m1n, m1z 158 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw 114 159 #endif 115 160 … … 148 193 ! 149 194 !-- Send left boundary, receive right boundary 150 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, &151 sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, &195 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, & 196 sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, & 152 197 comm2d, status, ierr ) 153 198 ! 154 199 !-- Send right boundary, receive left boundary 155 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, &156 sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, &200 CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, & 201 sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, & 157 202 comm2d, status, ierr ) 158 203 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) … … 192 237 DO j = nys, nyn 193 238 DO k = nzb+1, nzt 194 zaehler= ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )195 nenner= ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )196 fmax_l(1) = MAX( fmax_l(1) , zaehler)197 fmax_l(2) = MAX( fmax_l(2) , nenner)239 numera = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) ) 240 denomi = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 241 fmax_l(1) = MAX( fmax_l(1) , numera ) 242 fmax_l(2) = MAX( fmax_l(2) , denomi ) 198 243 ENDDO 199 244 ENDDO … … 210 255 ! 211 256 !-- Allocate temporary arrays 212 ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), &213 a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), &214 a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), &215 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &216 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &217 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &218 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), &219 sw(nzb+1:nzt,nxl-1:nxr+1) &257 ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), & 258 a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), & 259 a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), & 260 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), & 261 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), & 262 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), & 263 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), & 264 sw(nzb+1:nzt,nxl-1:nxr+1) & 220 265 ) 221 266 imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 … … 236 281 DO k = nzb+1, nzt 237 282 a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) 238 a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) &283 a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) & 239 284 + sk_p(k,j,i-1) ) 240 a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1) &241 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1) &285 a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1) & 286 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1) & 242 287 + 9.0 * sk_p(k,j,i-2) ) * f1920 243 a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1) &244 - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2) &288 a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1) & 289 - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2) & 245 290 ) * f48 246 a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1) &247 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1) &291 a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1) & 292 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1) & 248 293 - 3.0 * sk_p(k,j,i-2) ) * f48 249 294 ENDDO … … 259 304 cipf = 1.0 - 2.0 * cip 260 305 cimf = 1.0 - 2.0 * cim 261 ip = a0(k,i) * f2 * ( 1.0 - cipf ) &262 + a1(k,i) * f8 * ( 1.0 - cipf*cipf ) &306 ip = a0(k,i) * f2 * ( 1.0 - cipf ) & 307 + a1(k,i) * f8 * ( 1.0 - cipf*cipf ) & 263 308 + a2(k,i) * f24 * ( 1.0 - cipf*cipf*cipf ) 264 im = a0(k,i+1) * f2 * ( 1.0 - cimf ) &265 - a1(k,i+1) * f8 * ( 1.0 - cimf*cimf ) &309 im = a0(k,i+1) * f2 * ( 1.0 - cimf ) & 310 - a1(k,i+1) * f8 * ( 1.0 - cimf*cimf ) & 266 311 + a2(k,i+1) * f24 * ( 1.0 - cimf*cimf*cimf ) 267 312 ip = MAX( ip, 0.0 ) … … 274 319 cipf = 1.0 - 2.0 * cip 275 320 cimf = 1.0 - 2.0 * cim 276 ip = a0(k,i-1) * f2 * ( 1.0 - cipf ) &277 + a1(k,i-1) * f8 * ( 1.0 - cipf*cipf ) &321 ip = a0(k,i-1) * f2 * ( 1.0 - cipf ) & 322 + a1(k,i-1) * f8 * ( 1.0 - cipf*cipf ) & 278 323 + a2(k,i-1) * f24 * ( 1.0 - cipf*cipf*cipf ) 279 im = a0(k,i) * f2 * ( 1.0 - cimf ) &280 - a1(k,i) * f8 * ( 1.0 - cimf*cimf ) &324 im = a0(k,i) * f2 * ( 1.0 - cimf ) & 325 - a1(k,i) * f8 * ( 1.0 - cimf*cimf ) & 281 326 + a2(k,i) * f24 * ( 1.0 - cimf*cimf*cimf ) 282 327 ip = MAX( ip, 0.0 ) … … 309 354 DO i = nxl-1, nxr+1 310 355 DO k = nzb+1, nzt 311 m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / &356 m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / & 312 357 MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 ) 313 358 IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0 314 359 315 m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) / &360 m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) / & 316 361 MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35 ) 317 362 IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0 … … 322 367 323 368 !-- *VOCL STMT,IF(10) 324 IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0 &325 .OR. m2 > t2 .OR. m3 > T2 .OR.&326 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND. &327 m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 ) &369 IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0 & 370 .OR. m2 > t2 .OR. m3 > t2 .OR. & 371 ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND. & 372 m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 ) & 328 373 ) sw(k,i) = 1.0 329 374 ENDDO … … 425 470 DO i = nxl, nxr 426 471 DO k = nzb+1, nzt 427 fplus = ( 1.0 - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) &472 fplus = ( 1.0 - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) & 428 473 - ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i) 429 fminus = ( 1.0 - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) &474 fminus = ( 1.0 - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) & 430 475 - ( 1.0 - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i) 431 tend enz= fplus - fminus476 tendcy = fplus - fminus 432 477 ! 433 478 !-- Removed in order to optimize speed 434 479 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) 435 ! IF ( ( ABS( tend enz ) / ffmax ) < 1E-7 ) tendenz= 0.0480 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7 ) tendcy = 0.0 436 481 ! 437 482 !-- Density correction because of possible remaining divergences 438 483 d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx 439 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tend enz ) /&484 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 440 485 ( 1.0 + d_new ) 441 486 d(k,j,i) = d_new … … 447 492 ! 448 493 !-- Deallocate temporary arrays 449 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &494 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 450 495 ippb, ippe, m1, sw ) 451 496 … … 455 500 #if defined( __parallel ) 456 501 ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) 457 CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, &502 CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, & 458 503 type_xz_2, ierr ) 459 504 CALL MPI_TYPE_COMMIT( type_xz_2, ierr ) … … 461 506 !-- Send front boundary, receive rear boundary 462 507 CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) 463 CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, &464 sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, &508 CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, & 509 sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, & 465 510 comm2d, status, ierr ) 466 511 ! 467 512 !-- Send rear boundary, receive front boundary 468 CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, &469 sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, &513 CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, & 514 sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, & 470 515 comm2d, status, ierr ) 471 516 CALL MPI_TYPE_FREE( type_xz_2, ierr ) … … 490 535 DO j = nys, nyn 491 536 DO k = nzb+1, nzt 492 zaehler= ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )493 nenner= ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )494 fmax_l(1) = MAX( fmax_l(1) , zaehler)495 fmax_l(2) = MAX( fmax_l(2) , nenner)537 numera = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) ) 538 denomi = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 539 fmax_l(1) = MAX( fmax_l(1) , numera ) 540 fmax_l(2) = MAX( fmax_l(2) , denomi ) 496 541 ENDDO 497 542 ENDDO … … 508 553 ! 509 554 !-- Allocate temporary arrays 510 ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), &511 a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), &512 a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), &513 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), &514 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), &515 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), &516 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), &517 sw(nzb+1:nzt,nys-1:nyn+1) &555 ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), & 556 a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), & 557 a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), & 558 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), & 559 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), & 560 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), & 561 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), & 562 sw(nzb+1:nzt,nys-1:nyn+1) & 518 563 ) 519 564 imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 … … 528 573 DO k = nzb+1, nzt 529 574 a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) 530 a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) &575 a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) & 531 576 + sk_p(k,j-1,i) ) 532 a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i) &533 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i) &577 a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i) & 578 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i) & 534 579 + 9.0 * sk_p(k,j-2,i) ) * f1920 535 a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i) &536 - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i) &580 a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i) & 581 - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i) & 537 582 ) * f48 538 a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i) &539 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i) &583 a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i) & 584 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i) & 540 585 - 3.0 * sk_p(k,j-2,i) ) * f48 541 586 ENDDO … … 551 596 cipf = 1.0 - 2.0 * cip 552 597 cimf = 1.0 - 2.0 * cim 553 ip = a0(k,j) * f2 * ( 1.0 - cipf ) &554 + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &598 ip = a0(k,j) * f2 * ( 1.0 - cipf ) & 599 + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) & 555 600 + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf ) 556 im = a0(k,j+1) * f2 * ( 1.0 - cimf ) &557 - a1(k,j+1) * f8 * ( 1.0 - cimf*cimf ) &601 im = a0(k,j+1) * f2 * ( 1.0 - cimf ) & 602 - a1(k,j+1) * f8 * ( 1.0 - cimf*cimf ) & 558 603 + a2(k,j+1) * f24 * ( 1.0 - cimf*cimf*cimf ) 559 604 ip = MAX( ip, 0.0 ) … … 566 611 cipf = 1.0 - 2.0 * cip 567 612 cimf = 1.0 - 2.0 * cim 568 ip = a0(k,j-1) * f2 * ( 1.0 - cipf ) &569 + a1(k,j-1) * f8 * ( 1.0 - cipf*cipf ) &613 ip = a0(k,j-1) * f2 * ( 1.0 - cipf ) & 614 + a1(k,j-1) * f8 * ( 1.0 - cipf*cipf ) & 570 615 + a2(k,j-1) * f24 * ( 1.0 - cipf*cipf*cipf ) 571 im = a0(k,j) * f2 * ( 1.0 - cimf ) &572 - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &616 im = a0(k,j) * f2 * ( 1.0 - cimf ) & 617 - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) & 573 618 + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf ) 574 619 ip = MAX( ip, 0.0 ) … … 601 646 DO j = nys-1, nyn+1 602 647 DO k = nzb+1, nzt 603 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &648 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / & 604 649 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 ) 605 650 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 606 651 607 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &652 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / & 608 653 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 ) 609 654 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 … … 614 659 615 660 !-- *VOCL STMT,IF(10) 616 IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0 &617 .OR. m2 > t2 .OR. m3 > T2 .OR.&618 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND. &619 m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 ) &661 IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0 & 662 .OR. m2 > t2 .OR. m3 > t2 .OR. & 663 ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND. & 664 m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 ) & 620 665 ) sw(k,j) = 1.0 621 666 ENDDO … … 717 762 DO j = nys, nyn 718 763 DO k = nzb+1, nzt 719 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &764 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 720 765 - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j) 721 fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &766 fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) & 722 767 - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 723 tend enz= fplus - fminus768 tendcy = fplus - fminus 724 769 ! 725 770 !-- Removed in order to optimise speed 726 771 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) 727 ! IF ( ( ABS( tend enz ) / ffmax ) < 1E-7 ) tendenz= 0.0772 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7 ) tendcy = 0.0 728 773 ! 729 774 !-- Density correction because of possible remaining divergences 730 775 d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy 731 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tend enz ) /&776 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 732 777 ( 1.0 + d_new ) 733 778 d(k,j,i) = d_new … … 741 786 ! 742 787 !-- Deallocate temporary arrays 743 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &788 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 744 789 ippb, ippe, m1, sw ) 745 790 … … 879 924 ELSE 880 925 881 WRITE( message_string, * ) 'no vertical boundary condi', &926 WRITE( message_string, * ) 'no vertical boundary condi', & 882 927 'tion for variable "', sk_char, '"' 883 928 CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 ) … … 891 936 DO j = nys, nyn 892 937 DO k = nzb, nzt+1 893 zaehler= ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )894 nenner= ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )895 fmax_l(1) = MAX( fmax_l(1) , zaehler)896 fmax_l(2) = MAX( fmax_l(2) , nenner)938 numera = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) ) 939 denomi = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) ) 940 fmax_l(1) = MAX( fmax_l(1) , numera ) 941 fmax_l(2) = MAX( fmax_l(2) , denomi ) 897 942 ENDDO 898 943 ENDDO … … 909 954 ! 910 955 !-- Allocate temporary arrays 911 ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), &912 a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), &913 a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), &914 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), &915 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), &916 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), &917 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), &918 sw(nzb:nzt+1,nys:nyn) &956 ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), & 957 a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), & 958 a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), & 959 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), & 960 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), & 961 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), & 962 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), & 963 sw(nzb:nzt+1,nys:nyn) & 919 964 ) 920 965 imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 … … 929 974 DO k = nzb, nzt+1 930 975 a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) 931 a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &976 a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) & 932 977 + sk_p(k-1,j,i) ) 933 a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) &934 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) &978 a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) & 979 + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) & 935 980 + 9.0 * sk_p(k-2,j,i) ) * f1920 936 a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) &937 - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) &981 a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) & 982 - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) & 938 983 ) * f48 939 a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) &940 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) &984 a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) & 985 - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) & 941 986 - 3.0 * sk_p(k-2,j,i) ) * f48 942 987 ENDDO … … 952 997 cipf = 1.0 - 2.0 * cip 953 998 cimf = 1.0 - 2.0 * cim 954 ip = a0(k,j) * f2 * ( 1.0 - cipf ) &955 + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &999 ip = a0(k,j) * f2 * ( 1.0 - cipf ) & 1000 + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) & 956 1001 + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf ) 957 im = a0(k+1,j) * f2 * ( 1.0 - cimf ) &958 - a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) &1002 im = a0(k+1,j) * f2 * ( 1.0 - cimf ) & 1003 - a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) & 959 1004 + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf ) 960 1005 ip = MAX( ip, 0.0 ) … … 967 1012 cipf = 1.0 - 2.0 * cip 968 1013 cimf = 1.0 - 2.0 * cim 969 ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) &970 + a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) &1014 ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) & 1015 + a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) & 971 1016 + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf ) 972 im = a0(k,j) * f2 * ( 1.0 - cimf ) &973 - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &1017 im = a0(k,j) * f2 * ( 1.0 - cimf ) & 1018 - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) & 974 1019 + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf ) 975 1020 ip = MAX( ip, 0.0 ) … … 1002 1047 DO j = nys, nyn 1003 1048 DO k = nzb, nzt+1 1004 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &1049 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / & 1005 1050 MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 ) 1006 1051 IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 1007 1052 1008 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &1053 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / & 1009 1054 MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 ) 1010 1055 IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 … … 1015 1060 1016 1061 !-- *VOCL STMT,IF(10) 1017 IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 &1018 .OR. m2 > t2 .OR. m3 > T2 .OR.&1019 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. &1020 m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) &1062 IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 & 1063 .OR. m2 > t2 .OR. m3 > t2 .OR. & 1064 ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. & 1065 m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) & 1021 1066 ) sw(k,j) = 1.0 1022 1067 ENDDO … … 1118 1163 DO j = nys, nyn 1119 1164 DO k = nzb+1, nzt 1120 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &1165 fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & 1121 1166 - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j) 1122 fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &1167 fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) & 1123 1168 - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) 1124 tend enz= fplus - fminus1169 tendcy = fplus - fminus 1125 1170 ! 1126 1171 !-- Removed in order to optimise speed 1127 1172 ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) 1128 ! IF ( ( ABS( tend enz ) / ffmax ) < 1E-7 ) tendenz= 0.01173 ! IF ( ( ABS( tendcy ) / ffmax ) < 1E-7 ) tendcy = 0.0 1129 1174 ! 1130 1175 !-- Density correction because of possible remaining divergences 1131 1176 d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k) 1132 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tend enz ) /&1177 sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) / & 1133 1178 ( 1.0 + d_new ) 1134 1179 ! … … 1145 1190 DO j = nys, nyn 1146 1191 DO k = nzb+1, nzt 1147 sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &1192 sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + & 1148 1193 m1(k,j) * rmask(j,i,sr) 1149 1194 ENDDO … … 1158 1203 ! 1159 1204 !-- Deallocate temporary arrays 1160 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &1205 DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & 1161 1206 ippb, ippe, m1, sw ) 1162 1207
Note: See TracChangeset
for help on using the changeset viewer.
