SUBROUTINE advec_s_bc( sk, sk_char ) !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: advec_s_bc.f90 392 2009-09-24 10:39:14Z maronga $ ! ! 247 2009-02-27 14:01:30Z heinze ! Output of messages replaced by message handling routine ! ! 216 2008-11-25 07:12:43Z raasch ! Neumann boundary condition at k=nzb is explicitly set for better reading, ! although this has been already done in boundary_conds ! ! 97 2007-06-21 08:23:15Z raasch ! Advection of salinity included ! Bugfix: Error in boundary condition for TKE removed ! ! 63 2007-03-13 03:52:49Z raasch ! Calculation extended for gridpoint nzt ! ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.22 2006/02/23 09:42:08 raasch ! anz renamed ngp ! ! Revision 1.1 1997/08/29 08:53:46 raasch ! Initial revision ! ! ! Description: ! ------------ ! Advection term for scalar quantities using the Bott-Chlond scheme. ! Computation in individual steps for each of the three dimensions. ! Limiting assumptions: ! So far the scheme has been assuming equidistant grid spacing. As this is not ! the case in the stretched portion of the z-direction, there dzw(k) is used as ! a substitute for a constant grid length. This certainly causes incorrect ! results; however, it is hoped that they are not too apparent for weakly ! stretched grids. ! NOTE: This is a provisional, non-optimised version! !------------------------------------------------------------------------------! USE advection USE arrays_3d USE control_parameters USE cpulog USE grid_variables USE indices USE interfaces USE pegrid USE statistics IMPLICIT NONE CHARACTER (LEN=*) :: sk_char INTEGER :: i, ix, j, k, ngp, sr, type_xz_2 REAL :: cim, cimf, cip, cipf, d_new, ffmax, fminus, fplus, f2, f4, f8, & f12, f24, f48, f1920, im, ip, m2, m3, nenner, snenn, sterm, & tendenz, t1, t2, zaehler REAL :: fmax(2), fmax_l(2) REAL, DIMENSION(:,:,:), POINTER :: sk REAL, DIMENSION(:,:), ALLOCATABLE :: a0, a1, a12, a2, a22, immb, imme, & impb, impe, ipmb, ipme, ippb, ippe REAL, DIMENSION(:,:,:), ALLOCATABLE :: sk_p #if defined( __nec ) REAL (kind=4) :: m1n, m1z !Wichtig: Division REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw #else REAL :: m1n, m1z REAL, DIMENSION(:,:), ALLOCATABLE :: m1, sw #endif ! !-- Array sk_p requires 2 extra elements for each dimension ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) ) sk_p = 0.0 ! !-- Assign reciprocal values in order to avoid divisions later f2 = 0.5 f4 = 0.25 f8 = 0.125 f12 = 0.8333333333333333E-01 f24 = 0.4166666666666666E-01 f48 = 0.2083333333333333E-01 f1920 = 0.5208333333333333E-03 ! !-- Advection in x-direction: ! !-- Save the quantity to be advected in a local array !-- add an enlarged boundary in x-direction DO i = nxl-1, nxr+1 DO j = nys, nyn DO k = nzb, nzt+1 sk_p(k,j,i) = sk(k,j,i) ENDDO ENDDO ENDDO #if defined( __parallel ) ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' ) ! !-- Send left boundary, receive right boundary CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, & sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, & comm2d, status, ierr ) ! !-- Send right boundary, receive left boundary CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, & sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, & comm2d, status, ierr ) CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) #else ! !-- Cyclic boundary conditions sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2) sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1) sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1) sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2) #endif ! !-- In case of a sloping surface, the additional gridpoints in x-direction !-- of the temperature field at the left and right boundary of the total !-- domain must be adjusted by the temperature difference between this distance IF ( sloping_surface .AND. sk_char == 'pt' ) THEN IF ( nxl == 0 ) THEN sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset ENDIF IF ( nxr == nx ) THEN sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset ENDIF ENDIF ! !-- Initialise control density d = 0.0 ! !-- Determine maxima of the first and second derivative in x-direction fmax_l = 0.0 DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt zaehler = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) ) nenner = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) fmax_l(1) = MAX( fmax_l(1) , zaehler ) fmax_l(2) = MAX( fmax_l(2) , nenner ) ENDDO ENDDO ENDDO #if defined( __parallel ) CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) #else fmax = fmax_l #endif fmax = 0.04 * fmax ! !-- Allocate temporary arrays 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), a12(nzb+1:nzt,nxl-1:nxr+1), & a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), & imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), & impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), & ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), & ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), & sw(nzb+1:nzt,nxl-1:nxr+1) & ) imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 ! !-- Initialise point of time measuring of the exponential portion (this would !-- not work if done locally within the loop) CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' ) CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) ! !-- Outer loop of all j DO j = nys, nyn ! !-- Compute polynomial coefficients DO i = nxl-1, nxr+1 DO k = nzb+1, nzt a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) & + sk_p(k,j,i-1) ) a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1) & + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1) & + 9.0 * sk_p(k,j,i-2) ) * f1920 a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1) & - 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2) & ) * f48 a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1) & - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1) & - 3.0 * sk_p(k,j,i-2) ) * f48 ENDDO ENDDO ! !-- Fluxes using the Bott scheme !-- *VOCL LOOP,UNROLL(2) DO i = nxl, nxr DO k = nzb+1, nzt cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k,i) * f2 * ( 1.0 - cipf ) & + a1(k,i) * f8 * ( 1.0 - cipf*cipf ) & + a2(k,i) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k,i+1) * f2 * ( 1.0 - cimf ) & - a1(k,i+1) * f8 * ( 1.0 - cimf*cimf ) & + a2(k,i+1) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ippb(k,i) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) impb(k,i) = im * MIN( 1.0, sk_p(k,j,i+1) / (ip+im+1E-15) ) cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k,i-1) * f2 * ( 1.0 - cipf ) & + a1(k,i-1) * f8 * ( 1.0 - cipf*cipf ) & + a2(k,i-1) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k,i) * f2 * ( 1.0 - cimf ) & - a1(k,i) * f8 * ( 1.0 - cimf*cimf ) & + a2(k,i) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ipmb(k,i) = ip * MIN( 1.0, sk_p(k,j,i-1) / (ip+im+1E-15) ) immb(k,i) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) ENDDO ENDDO ! !-- Compute monitor function m1 DO i = nxl-2, nxr+2 DO k = nzb+1, nzt m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) ) m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) ) IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN m1(k,i) = m1z / m1n IF ( m1(k,i) /= 2.0 .AND. m1n < fmax(2) ) m1(k,i) = 0.0 ELSEIF ( m1n < m1z ) THEN m1(k,i) = -1.0 ELSE m1(k,i) = 0.0 ENDIF ENDDO ENDDO ! !-- Compute switch sw sw = 0.0 DO i = nxl-1, nxr+1 DO k = nzb+1, nzt m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / & MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 ) IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0 m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) / & MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35 ) IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0 t1 = 0.35 t2 = 0.35 IF ( m1(k,i) == -1.0 ) t2 = 0.12 !-- *VOCL STMT,IF(10) IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0 & .OR. m2 > t2 .OR. m3 > T2 .OR. & ( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND. & m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 ) & ) sw(k,i) = 1.0 ENDDO ENDDO ! !-- Fluxes using the exponential scheme CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) DO i = nxl, nxr DO k = nzb+1, nzt !-- *VOCL STMT,IF(10) IF ( sw(k,i) == 1.0 ) THEN snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ippe(k,i) = sk_p(k,j,i) * cip IF ( sterm == 0.9999 ) ippe(k,i) = sk_p(k,j,i) * cip snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) imme(k,i) = sk_p(k,j,i) * cim IF ( sterm == 0.9999 ) imme(k,i) = sk_p(k,j,i) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k,i+1) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k,j,i+2) IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx ) impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) impe(k,i) = sk_p(k,j,i+1) * cim IF ( sterm == 0.9999 ) impe(k,i) = sk_p(k,j,i+1) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k,i-1) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k,j,i-2) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx ) ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ipme(k,i) = sk_p(k,j,i-1) * cip IF ( sterm == 0.9999 ) ipme(k,i) = sk_p(k,j,i-1) * cip ENDIF ENDDO ENDDO CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) ! !-- Prognostic equation DO i = nxl, nxr DO k = nzb+1, nzt fplus = ( 1.0 - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) & - ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i) fminus = ( 1.0 - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) & - ( 1.0 - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i) tendenz = fplus - fminus ! !-- Removed in order to optimize speed ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) ! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0 ! !-- Density correction because of possible remaining divergences d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / & ( 1.0 + d_new ) d(k,j,i) = d_new ENDDO ENDDO ENDDO ! End of the advection in x-direction ! !-- Deallocate temporary arrays DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & ippb, ippe, m1, sw ) ! !-- Enlarge boundary of local array cyclically in y-direction #if defined( __parallel ) ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 ) CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, & type_xz_2, ierr ) CALL MPI_TYPE_COMMIT( type_xz_2, ierr ) ! !-- Send front boundary, receive rear boundary CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, & sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, & comm2d, status, ierr ) ! !-- Send rear boundary, receive front boundary CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, & sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, & comm2d, status, ierr ) CALL MPI_TYPE_FREE( type_xz_2, ierr ) CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' ) #else DO i = nxl, nxr DO k = nzb+1, nzt sk_p(k,nys-1,i) = sk_p(k,nyn,i) sk_p(k,nys-2,i) = sk_p(k,nyn-1,i) sk_p(k,nys-3,i) = sk_p(k,nyn-2,i) sk_p(k,nyn+1,i) = sk_p(k,nys,i) sk_p(k,nyn+2,i) = sk_p(k,nys+1,i) sk_p(k,nyn+3,i) = sk_p(k,nys+2,i) ENDDO ENDDO #endif ! !-- Determine the maxima of the first and second derivative in y-direction fmax_l = 0.0 DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt zaehler = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) ) nenner = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) fmax_l(1) = MAX( fmax_l(1) , zaehler ) fmax_l(2) = MAX( fmax_l(2) , nenner ) ENDDO ENDDO ENDDO #if defined( __parallel ) CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) #else fmax = fmax_l #endif fmax = 0.04 * fmax ! !-- Allocate temporary arrays 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), a12(nzb+1:nzt,nys-1:nyn+1), & a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), & imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), & impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), & ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), & ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), & sw(nzb+1:nzt,nys-1:nyn+1) & ) imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 ! !-- Outer loop of all i DO i = nxl, nxr ! !-- Compute polynomial coefficients DO j = nys-1, nyn+1 DO k = nzb+1, nzt a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) & + sk_p(k,j-1,i) ) a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i) & + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i) & + 9.0 * sk_p(k,j-2,i) ) * f1920 a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i) & - 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i) & ) * f48 a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i) & - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i) & - 3.0 * sk_p(k,j-2,i) ) * f48 ENDDO ENDDO ! !-- Fluxes using the Bott scheme !-- *VOCL LOOP,UNROLL(2) DO j = nys, nyn DO k = nzb+1, nzt cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k,j) * f2 * ( 1.0 - cipf ) & + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) & + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k,j+1) * f2 * ( 1.0 - cimf ) & - a1(k,j+1) * f8 * ( 1.0 - cimf*cimf ) & + a2(k,j+1) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) impb(k,j) = im * MIN( 1.0, sk_p(k,j+1,i) / (ip+im+1E-15) ) cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k,j-1) * f2 * ( 1.0 - cipf ) & + a1(k,j-1) * f8 * ( 1.0 - cipf*cipf ) & + a2(k,j-1) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k,j) * f2 * ( 1.0 - cimf ) & - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) & + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ipmb(k,j) = ip * MIN( 1.0, sk_p(k,j-1,i) / (ip+im+1E-15) ) immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) ENDDO ENDDO ! !-- Compute monitor function m1 DO j = nys-2, nyn+2 DO k = nzb+1, nzt m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) ) m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) ) IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN m1(k,j) = m1z / m1n IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0 ELSEIF ( m1n < m1z ) THEN m1(k,j) = -1.0 ELSE m1(k,j) = 0.0 ENDIF ENDDO ENDDO ! !-- Compute switch sw sw = 0.0 DO j = nys-1, nyn+1 DO k = nzb+1, nzt m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / & MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 ) IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / & MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 ) IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 t1 = 0.35 t2 = 0.35 IF ( m1(k,j) == -1.0 ) t2 = 0.12 !-- *VOCL STMT,IF(10) IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0 & .OR. m2 > t2 .OR. m3 > T2 .OR. & ( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND. & m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 ) & ) sw(k,j) = 1.0 ENDDO ENDDO ! !-- Fluxes using exponential scheme CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) DO j = nys, nyn DO k = nzb+1, nzt !-- *VOCL STMT,IF(10) IF ( sw(k,j) == 1.0 ) THEN snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k,j+1) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k,j+2,i) IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy ) impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k,j+1,i) * cim IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k,j+1,i) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k,j-1) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k,j-2,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy ) ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k,j-1,i) * cip IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k,j-1,i) * cip ENDIF ENDDO ENDDO CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) ! !-- Prognostic equation DO j = nys, nyn DO k = nzb+1, nzt fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j) fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) & - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) tendenz = fplus - fminus ! !-- Removed in order to optimise speed ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) ! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0 ! !-- Density correction because of possible remaining divergences d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / & ( 1.0 + d_new ) d(k,j,i) = d_new ENDDO ENDDO ENDDO ! End of the advection in y-direction CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' ) CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' ) ! !-- Deallocate temporary arrays DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & ippb, ippe, m1, sw ) ! !-- Initialise for the computation of heat fluxes (see below; required in !-- UP flow_statistics) IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0 ! !-- Add top and bottom boundaries according to the relevant boundary conditions IF ( sk_char == 'pt' ) THEN ! !-- Temperature boundary condition at the bottom boundary IF ( ibc_pt_b == 0 ) THEN ! !-- Dirichlet (fixed surface temperature) DO i = nxl, nxr DO j = nys, nyn sk_p(nzb-1,j,i) = sk_p(nzb,j,i) sk_p(nzb-2,j,i) = sk_p(nzb,j,i) ENDDO ENDDO ELSE ! !-- Neumann (i.e. here zero gradient) DO i = nxl, nxr DO j = nys, nyn sk_p(nzb,j,i) = sk_p(nzb+1,j,i) sk_p(nzb-1,j,i) = sk_p(nzb,j,i) sk_p(nzb-2,j,i) = sk_p(nzb,j,i) ENDDO ENDDO ENDIF ! !-- Temperature boundary condition at the top boundary IF ( ibc_pt_t == 0 .OR. ibc_pt_t == 1 ) THEN ! !-- Dirichlet or Neumann (zero gradient) DO i = nxl, nxr DO j = nys, nyn sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) ENDDO ENDDO ELSEIF ( ibc_pt_t == 2 ) THEN ! !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead DO i = nxl, nxr DO j = nys, nyn sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1) sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1) ENDDO ENDDO ENDIF ELSEIF ( sk_char == 'sa' ) THEN ! !-- Salinity boundary condition at the bottom boundary. !-- So far, always Neumann (i.e. here zero gradient) is used DO i = nxl, nxr DO j = nys, nyn sk_p(nzb,j,i) = sk_p(nzb+1,j,i) sk_p(nzb-1,j,i) = sk_p(nzb,j,i) sk_p(nzb-2,j,i) = sk_p(nzb,j,i) ENDDO ENDDO ! !-- Salinity boundary condition at the top boundary. !-- Dirichlet or Neumann (zero gradient) DO i = nxl, nxr DO j = nys, nyn sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) ENDDO ENDDO ELSEIF ( sk_char == 'q' ) THEN ! !-- Specific humidity boundary condition at the bottom boundary. !-- Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient) DO i = nxl, nxr DO j = nys, nyn sk_p(nzb,j,i) = sk_p(nzb+1,j,i) sk_p(nzb-1,j,i) = sk_p(nzb,j,i) sk_p(nzb-2,j,i) = sk_p(nzb,j,i) ENDDO ENDDO ! !-- Specific humidity boundary condition at the top boundary IF ( ibc_q_t == 0 ) THEN ! !-- Dirichlet DO i = nxl, nxr DO j = nys, nyn sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) ENDDO ENDDO ELSE ! !-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead DO i = nxl, nxr DO j = nys, nyn sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1) sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1) ENDDO ENDDO ENDIF ELSEIF ( sk_char == 'e' ) THEN ! !-- TKE boundary condition at bottom and top boundary (generally Neumann) DO i = nxl, nxr DO j = nys, nyn sk_p(nzb,j,i) = sk_p(nzb+1,j,i) sk_p(nzb-1,j,i) = sk_p(nzb,j,i) sk_p(nzb-2,j,i) = sk_p(nzb,j,i) sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) ENDDO ENDDO ELSE WRITE( message_string, * ) 'no vertical boundary condi', & 'tion for variable "', sk_char, '"' CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 ) ENDIF ! !-- Determine the maxima of the first and second derivative in z-direction fmax_l = 0.0 DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) ) nenner = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) ) fmax_l(1) = MAX( fmax_l(1) , zaehler ) fmax_l(2) = MAX( fmax_l(2) , nenner ) ENDDO ENDDO ENDDO #if defined( __parallel ) CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr ) #else fmax = fmax_l #endif fmax = 0.04 * fmax ! !-- Allocate temporary arrays ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), & a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), & a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), & imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), & impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), & ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), & ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), & sw(nzb:nzt+1,nys:nyn) & ) imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0 ! !-- Outer loop of all i DO i = nxl, nxr ! !-- Compute polynomial coefficients DO j = nys, nyn DO k = nzb, nzt+1 a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) & + sk_p(k-1,j,i) ) a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) & + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) & + 9.0 * sk_p(k-2,j,i) ) * f1920 a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) & - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) & ) * f48 a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) & - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) & - 3.0 * sk_p(k-2,j,i) ) * f48 ENDDO ENDDO ! !-- Fluxes using the Bott scheme !-- *VOCL LOOP,UNROLL(2) DO j = nys, nyn DO k = nzb+1, nzt cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) ) cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k,j) * f2 * ( 1.0 - cipf ) & + a1(k,j) * f8 * ( 1.0 - cipf*cipf ) & + a2(k,j) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k+1,j) * f2 * ( 1.0 - cimf ) & - a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) & + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) impb(k,j) = im * MIN( 1.0, sk_p(k+1,j,i) / (ip+im+1E-15) ) cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) ) cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) ) cipf = 1.0 - 2.0 * cip cimf = 1.0 - 2.0 * cim ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) & + a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) & + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf ) im = a0(k,j) * f2 * ( 1.0 - cimf ) & - a1(k,j) * f8 * ( 1.0 - cimf*cimf ) & + a2(k,j) * f24 * ( 1.0 - cimf*cimf*cimf ) ip = MAX( ip, 0.0 ) im = MAX( im, 0.0 ) ipmb(k,j) = ip * MIN( 1.0, sk_p(k-1,j,i) / (ip+im+1E-15) ) immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) ) ENDDO ENDDO ! !-- Compute monitor function m1 DO j = nys, nyn DO k = nzb-1, nzt+2 m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) ) m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) ) IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN m1(k,j) = m1z / m1n IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0 ELSEIF ( m1n < m1z ) THEN m1(k,j) = -1.0 ELSE m1(k,j) = 0.0 ENDIF ENDDO ENDDO ! !-- Compute switch sw sw = 0.0 DO j = nys, nyn DO k = nzb, nzt+1 m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / & MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 ) IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0 m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / & MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 ) IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0 t1 = 0.35 t2 = 0.35 IF ( m1(k,j) == -1.0 ) t2 = 0.12 !-- *VOCL STMT,IF(10) IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 & .OR. m2 > t2 .OR. m3 > T2 .OR. & ( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. & m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) & ) sw(k,j) = 1.0 ENDDO ENDDO ! !-- Fluxes using exponential scheme CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) DO j = nys, nyn DO k = nzb+1, nzt !-- *VOCL STMT,IF(10) IF ( sw(k,j) == 1.0 ) THEN snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) ) ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) ) imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k+1,j) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k+2,j,i) IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) ) impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( & aex(ix) * cim + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) & ) & ) IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k+1,j,i) * cim IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k+1,j,i) * cim ENDIF !-- *VOCL STMT,IF(10) IF ( sw(k-1,j) == 1.0 ) THEN snenn = sk_p(k,j,i) - sk_p(k-2,j,i) IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9 sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn sterm = MIN( sterm, 0.9999 ) sterm = MAX( sterm, 0.0001 ) ix = INT( sterm * 1000 ) + 1 cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) ) ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( & aex(ix) * cip + bex(ix) / dex(ix) * ( & eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) & ) & ) IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k-1,j,i) * cip IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k-1,j,i) * cip ENDIF ENDDO ENDDO CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' ) ! !-- Prognostic equation DO j = nys, nyn DO k = nzb+1, nzt fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) & - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j) fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) & - ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j) tendenz = fplus - fminus ! !-- Removed in order to optimise speed ! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 ) ! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0 ! !-- Density correction because of possible remaining divergences d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k) sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / & ( 1.0 + d_new ) ! !-- Store heat flux for subsequent statistics output. !-- array m1 is here used as temporary storage m1(k,j) = fplus / dt_3d * dzw(k) ENDDO ENDDO ! !-- Sum up heat flux in order to order to obtain horizontal averages IF ( sk_char == 'pt' ) THEN DO sr = 0, statistic_regions DO j = nys, nyn DO k = nzb+1, nzt sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + & m1(k,j) * rmask(j,i,sr) ENDDO ENDDO ENDDO ENDIF ENDDO ! End of the advection in z-direction CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' ) CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' ) ! !-- Deallocate temporary arrays DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, & ippb, ippe, m1, sw ) ! !-- Store results as tendency and deallocate local array DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d ENDDO ENDDO ENDDO DEALLOCATE( sk_p ) END SUBROUTINE advec_s_bc