- Timestamp:
- Nov 2, 2018 2:41:25 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
palm.f90 (modified) (2 diffs)
-
pmc_interface_mod.f90 (modified) (147 diffs)
-
time_integration.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/palm.f90
r3458 r3484 25 25 ! ----------------- 26 26 ! $Id$ 27 ! pmci_ensure_nest_mass_conservation removed permanently 28 ! 29 ! 3458 2018-10-30 14:51:23Z kanani 27 30 ! from chemistry branch r3443, forkel: 28 31 ! removed double do_emis check around CALL chem_init … … 308 311 USE pmc_interface, & 309 312 ONLY: nested_run, pmci_child_initialize, pmci_init, & 310 pmci_modelconfiguration, pmci_parent_initialize, & 311 pmci_ensure_nest_mass_conservation 313 pmci_modelconfiguration, pmci_parent_initialize 312 314 313 315 USE write_restart_data_mod, & -
palm/trunk/SOURCE/pmc_interface_mod.f90
r3274 r3484 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Introduction of reversibility correction to the interpolation routines in order to 28 ! guarantee mass and scalar conservation through the nest boundaries. Several errors 29 ! are corrected and pmci_ensure_nest_mass_conservation is permanently removed. 30 ! 31 ! 3274 2018-09-24 15:42:55Z knoop 27 32 ! Modularization of all bulk cloud physics code components 28 33 ! … … 366 371 USE pegrid, & 367 372 ONLY: collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy, & 368 numprocs 373 numprocs, pleft, pnorth, pright, psouth, status 369 374 370 375 USE pmc_child, & … … 417 422 INTEGER(iwp), SAVE :: cpl_parent_id !< 418 423 ! 419 !-- Control parameters , will be made input parameters later424 !-- Control parameters 420 425 CHARACTER(LEN=7), SAVE :: nesting_datatransfer_mode = 'mixed' !< steering 421 426 !< parameter for data- … … 441 446 442 447 ! 443 !-- Child coarse data arrays 444 INTEGER(iwp), DIMENSION(5),PUBLIC :: coarse_bound !< 448 !-- Children's parent-grid arrays 449 INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC :: coarse_bound !< subdomain index bounds for children's parent-grid arrays 450 INTEGER(iwp), SAVE, DIMENSION(4), PUBLIC :: coarse_bound_anterp !< subdomain index bounds for anterpolation 445 451 446 452 REAL(wp), SAVE :: xexl !< … … 492 498 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: r1zw !< 493 499 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: r2zw !< 500 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) :: celltmpd !< 494 501 ! 495 502 !-- Child index arrays and log-ratio arrays for the log-law near-wall … … 534 541 ! 535 542 !-- Upper bounds for k in anterpolation. 536 INTEGER(iwp), SAVE :: kct u!<543 INTEGER(iwp), SAVE :: kcto !< 537 544 INTEGER(iwp), SAVE :: kctw !< 538 545 ! … … 599 606 600 607 TYPE(coarsegrid_def), SAVE, PUBLIC :: cg !< 601 602 !- Variables for particle coupling 603 608 ! 609 !-- Variables for particle coupling 604 610 TYPE, PUBLIC :: childgrid_def 605 611 INTEGER(iwp) :: nx !< … … 641 647 END INTERFACE pmci_datatrans 642 648 643 INTERFACE pmci_ensure_nest_mass_conservation644 MODULE PROCEDURE pmci_ensure_nest_mass_conservation645 END INTERFACE646 647 649 INTERFACE pmci_init 648 650 MODULE PROCEDURE pmci_init … … 687 689 PUBLIC pmci_child_initialize 688 690 PUBLIC pmci_datatrans 689 PUBLIC pmci_ensure_nest_mass_conservation690 691 PUBLIC pmci_init 691 692 PUBLIC pmci_modelconfiguration … … 891 892 IF ( myid == 0 ) THEN 892 893 893 CALL pmc_recv_from_child( child_id, val, size(val), 0, 123, ierr )894 CALL pmc_recv_from_child( child_id, fval, size(fval), 0, 124, ierr )894 CALL pmc_recv_from_child( child_id, val, SIZE(val), 0, 123, ierr ) 895 CALL pmc_recv_from_child( child_id, fval, SIZE(fval), 0, 124, ierr ) 895 896 896 897 nx_cl = val(1) … … 983 984 mid = pmc_parent_for_child(mm) 984 985 ! 985 !-- Check Only different nest level986 !-- Check only different nest levels 986 987 IF (m_couplers(child_id)%parent_id /= m_couplers(mid)%parent_id) THEN 987 988 IF ( ( ch_xl(m) < ch_xr(mm) .OR. & … … 1000 1001 DEALLOCATE( cl_coord_x ) 1001 1002 DEALLOCATE( cl_coord_y ) 1002 1003 1003 ! 1004 1004 !-- Send information about operating mode (LES or RANS) to child. This will be … … 1219 1219 1220 1220 INTEGER(iwp) :: i !< 1221 INTEGER(iwp) :: ierr !< 1222 INTEGER(iwp) :: icl !< 1223 INTEGER(iwp) :: icr !< 1224 INTEGER(iwp) :: j !< 1225 INTEGER(iwp) :: jcn !< 1226 INTEGER(iwp) :: jcs !< 1221 INTEGER(iwp) :: ierr !< 1222 INTEGER(iwp) :: icl !< left index limit for children's parent-grid arrays 1223 INTEGER(iwp) :: icla !< left index limit for anterpolation 1224 INTEGER(iwp) :: icr !< left index limit for children's parent-grid arrays 1225 INTEGER(iwp) :: icra !< right index limit for anterpolation 1226 INTEGER(iwp) :: j !< 1227 INTEGER(iwp) :: jcn !< north index limit for children's parent-grid arrays 1228 INTEGER(iwp) :: jcna !< north index limit for anterpolation 1229 INTEGER(iwp) :: jcs !< sout index limit for children's parent-grid arrays 1230 INTEGER(iwp) :: jcsa !< south index limit for anterpolation 1227 1231 INTEGER(iwp) :: n !< running index for number of chemical species 1228 1232 … … 1374 1378 !-- Get coarse grid coordinates and values of the z-direction from the parent 1375 1379 IF ( myid == 0) THEN 1376 1377 1380 CALL pmc_recv_from_parent( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr ) 1378 1381 CALL pmc_recv_from_parent( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr ) 1379 CALL pmc_recv_from_parent( cg%dzu, cg%nz + 1, 0, 26, ierr ) 1380 CALL pmc_recv_from_parent( cg%dzw, cg%nz + 1, 0, 27, ierr ) 1381 CALL pmc_recv_from_parent( cg%zu, cg%nz + 2, 0, 28, ierr ) 1382 CALL pmc_recv_from_parent( cg%zw, cg%nz + 2, 0, 29, ierr ) 1383 1382 CALL pmc_recv_from_parent( cg%dzu, cg%nz+1, 0, 26, ierr ) 1383 CALL pmc_recv_from_parent( cg%dzw, cg%nz+1, 0, 27, ierr ) 1384 CALL pmc_recv_from_parent( cg%zu, cg%nz+2, 0, 28, ierr ) 1385 CALL pmc_recv_from_parent( cg%zw, cg%nz+2, 0, 29, ierr ) 1384 1386 ENDIF 1385 1387 ! … … 1392 1394 CALL MPI_BCAST( cg%zw, cg%nz+2, MPI_REAL, 0, comm2d, ierr ) 1393 1395 CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr ) 1394 1395 1396 ! 1396 1397 !-- Find the index bounds for the nest domain in the coarse-grid index space … … 1441 1442 !-- Precompute the index arrays and relaxation functions for the 1442 1443 !-- anterpolation 1443 IF ( TRIM( nesting_mode ) == 'two-way' .OR. & 1444 nesting_mode == 'vertical' ) THEN 1445 CALL pmci_init_anterp_tophat 1446 ENDIF 1444 ! 1445 !-- Note that the anterpolation index bounds are needed also in case 1446 !-- of one-way coupling because of the reversibility correction 1447 !-- included in the interpolation algorithms. 1448 CALL pmci_init_anterp_tophat 1447 1449 ! 1448 1450 !-- Finally, compute the total area of the top-boundary face of the domain. … … 1463 1465 INTEGER(iwp), DIMENSION(5,numprocs) :: coarse_bound_all !< 1464 1466 INTEGER(iwp), DIMENSION(2) :: size_of_array !< 1465 1467 1468 INTEGER(iwp) :: i !< 1469 INTEGER(iwp) :: ijaux !< 1470 INTEGER(iwp) :: j !< 1466 1471 REAL(wp) :: loffset !< 1467 1472 REAL(wp) :: noffset !< … … 1536 1541 coarse_bound(5) = myid 1537 1542 ! 1543 !-- Determine the anterpolation index limits. If at least half of the 1544 !-- parent-grid cell is within the current child sub-domain, then it 1545 !-- is included in the current sub-domain's anterpolation domain. 1546 !-- Else the parent-grid cell is included in the neighbouring subdomain's 1547 !-- anterpolation domain, or not included at all if we are at the outer 1548 !-- edge of the child domain. 1549 DO i = 0, cg%nx 1550 IF ( cg%coord_x(i) + 0.5_wp * cg%dx >= coord_x(nxl) ) THEN 1551 icla = MAX( 0, i ) 1552 EXIT 1553 ENDIF 1554 ENDDO 1555 DO i = cg%nx, 0 , -1 1556 IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= coord_x(nxr+1) ) THEN 1557 icra = MIN( cg%nx, i ) 1558 EXIT 1559 ENDIF 1560 ENDDO 1561 DO j = 0, cg%ny 1562 IF ( cg%coord_y(j) + 0.5_wp * cg%dy >= coord_y(nys) ) THEN 1563 jcsa = MAX( 0, j ) 1564 EXIT 1565 ENDIF 1566 ENDDO 1567 DO j = cg%ny, 0 , -1 1568 IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= coord_y(nyn+1) ) THEN 1569 jcna = MIN( cg%ny, j ) 1570 EXIT 1571 ENDIF 1572 ENDDO 1573 ! 1574 !-- Make sure that the indexing is contiguous 1575 IF ( nxl == 0 ) THEN 1576 CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr ) 1577 ELSE IF ( nxr == nx ) THEN 1578 CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 1579 icla = ijaux + 1 1580 ELSE 1581 CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr ) 1582 CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 1583 icla = ijaux + 1 1584 ENDIF 1585 IF ( nys == 0 ) THEN 1586 CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr ) 1587 ELSE IF ( nyn == ny ) THEN 1588 CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 1589 jcsa = ijaux + 1 1590 ELSE 1591 CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr ) 1592 CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 1593 jcsa = ijaux + 1 1594 ENDIF 1595 1596 write(9,"('Anterpolation bounds: ',4(i3,2x))") icla, icra, jcsa, jcna 1597 flush(9) 1598 coarse_bound_anterp(1) = icla 1599 coarse_bound_anterp(2) = icra 1600 coarse_bound_anterp(3) = jcsa 1601 coarse_bound_anterp(4) = jcna 1602 ! 1538 1603 !-- Note that MPI_Gather receives data from all processes in the rank order 1539 1604 !-- TO_DO: refer to the line where this fact becomes important … … 1561 1626 IMPLICIT NONE 1562 1627 1628 INTEGER(iwp) :: acsize !< 1563 1629 INTEGER(iwp) :: i !< 1564 1630 INTEGER(iwp) :: j !< … … 1568 1634 INTEGER(iwp) :: kdzw !< 1569 1635 1636 REAL(wp) :: dzmin !< 1637 REAL(wp) :: parentdzmax !< 1570 1638 REAL(wp) :: xb !< 1571 1639 REAL(wp) :: xcsu !< … … 1583 1651 REAL(wp) :: zfsw !< 1584 1652 1585 1653 1586 1654 xb = nxl * dx 1587 1655 yb = nys * dy … … 1658 1726 r1zo(k) = 1.0_wp - r2zo(k) 1659 1727 ENDDO 1728 ! 1729 !-- Determine the maximum dimension of anterpolation cells and allocate the 1730 !-- work array celltmpd needed in the reversibility correction in the 1731 !-- interpolation 1732 dzmin = 999999.9_wp 1733 DO k = 1, nzt+1 1734 dzmin = MIN( dzmin, dzu(k), dzw(k) ) 1735 ENDDO 1736 parentdzmax = 0.0_wp 1737 DO k = 1, cg%nz+1 1738 parentdzmax = MAX(parentdzmax , cg%dzu(k), cg%dzw(k) ) 1739 ENDDO 1740 acsize = CEILING( cg%dx / dx ) * CEILING( cg%dy / dy ) * & 1741 CEILING( parentdzmax / dzmin ) 1742 ALLOCATE( celltmpd(1:acsize) ) 1743 ! write(9,"('acsize: ',i3,2(e12.5,2x))") acsize, dzmin, parentdzmax 1660 1744 1661 1745 END SUBROUTINE pmci_init_interp_tril … … 3007 3091 REAL(wp) :: xi !< 3008 3092 REAL(wp) :: eta !< 3093 REAL(wp) :: tolerance !< 3009 3094 REAL(wp) :: zeta !< 3010 3095 … … 3027 3112 ENDIF 3028 3113 ! 3029 !-- First determine kct uand kctw that are the coarse-grid upper bounds for3114 !-- First determine kcto and kctw that are the coarse-grid upper bounds for 3030 3115 !-- index k 3031 3116 kk = 0 … … 3033 3118 kk = kk + 1 3034 3119 ENDDO 3035 kct u= kk - 13120 kcto = kk - 1 3036 3121 3037 3122 kk = 0 … … 3049 3134 ALLOCATE( jfuv(jcs:jcn) ) 3050 3135 ALLOCATE( jfuo(jcs:jcn) ) 3051 ALLOCATE( kflw(0:kctw) ) 3052 ALLOCATE( kflo(0:kctu) ) 3053 ALLOCATE( kfuw(0:kctw) ) 3054 ALLOCATE( kfuo(0:kctu) ) 3055 3056 ALLOCATE( ijkfc_u(0:kctu,jcs:jcn,icl:icr) ) 3057 ALLOCATE( ijkfc_v(0:kctu,jcs:jcn,icl:icr) ) 3058 ALLOCATE( ijkfc_w(0:kctw,jcs:jcn,icl:icr) ) 3059 ALLOCATE( ijkfc_s(0:kctu,jcs:jcn,icl:icr) ) 3136 ALLOCATE( kflw(0:cg%nz+1) ) 3137 ALLOCATE( kflo(0:cg%nz+1) ) 3138 ALLOCATE( kfuw(0:cg%nz+1) ) 3139 ALLOCATE( kfuo(0:cg%nz+1) ) 3140 3141 ALLOCATE( ijkfc_u(0:cg%nz+1,jcs:jcn,icl:icr) ) 3142 ALLOCATE( ijkfc_v(0:cg%nz+1,jcs:jcn,icl:icr) ) 3143 ALLOCATE( ijkfc_w(0:cg%nz+1,jcs:jcn,icl:icr) ) 3144 ALLOCATE( ijkfc_s(0:cg%nz+1,jcs:jcn,icl:icr) ) 3145 3060 3146 ijkfc_u = 0 3061 3147 ijkfc_v = 0 … … 3065 3151 !-- i-indices of u for each ii-index value 3066 3152 !-- ii=icr is redundant for anterpolation 3153 tolerance = 0.000001_wp * dx 3067 3154 istart = nxlg 3068 3155 DO ii = icl, icr-1 3156 ! 3157 !-- In case the child and parent grid lines match in x 3158 !-- use only the local k,j-child-grid plane for the anterpolation, 3159 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3069 3160 i = istart 3070 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) .AND. & 3071 ( i < nxrg ) ) 3072 i = i + 1 3161 DO WHILE ( coord_x(i) < cg%coord_x(ii) - tolerance .AND. i < nxrg ) 3162 i = i + 1 3073 3163 ENDDO 3074 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3075 ir = i 3076 DO WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx ) .AND.& 3077 ( i < nxrg+1 ) ) 3078 i = i + 1 3079 ir = MIN( i, nxrg ) 3080 ENDDO 3081 ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg ) 3082 istart = iflu(ii) 3164 IF ( ABS( coord_x(i) - cg%coord_x(ii) ) < tolerance ) THEN 3165 i = istart 3166 DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg ) 3167 i = i + 1 3168 ENDDO 3169 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3170 ifuu(ii) = iflu(ii) 3171 istart = iflu(ii) 3172 ELSE 3173 i = istart 3174 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) & 3175 .AND. ( i < nxrg ) ) 3176 i = i + 1 3177 ENDDO 3178 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3179 ir = i 3180 DO WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx ) & 3181 .AND. ( i < nxrg+1 ) ) 3182 i = i + 1 3183 ir = MIN( i, nxrg ) 3184 ENDDO 3185 ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg ) 3186 istart = iflu(ii) 3187 ENDIF 3188 !AH 3189 write(9,"('pmci_init_anterp_tophat, ii, iflu, ifuu: ', 3(i4,2x))") & 3190 ii, iflu(ii), ifuu(ii) 3191 flush(9) 3192 3083 3193 ENDDO 3084 3194 iflu(icr) = nxrg … … 3104 3214 istart = iflo(ii) 3105 3215 ENDDO 3216 !AH 3217 write(9,"('pmci_init_anterp_tophat, ii, iflo, ifuo: ', 3(i4,2x))") & 3218 ii, iflo(ii), ifuo(ii) 3219 flush(9) 3220 3106 3221 iflo(icr) = nxrg 3107 3222 ifuo(icr) = nxrg … … 3109 3224 !-- j-indices of v for each jj-index value 3110 3225 !-- jj=jcn is redundant for anterpolation 3226 tolerance = 0.000001_wp * dy 3111 3227 jstart = nysg 3112 3228 DO jj = jcs, jcn-1 3229 ! 3230 !-- In case the child and parent grid lines match in y 3231 !-- use only the local k,i-child-grid plane for the anterpolation, 3232 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3113 3233 j = jstart 3114 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) .AND. & 3115 ( j < nyng ) ) 3116 j = j + 1 3234 DO WHILE ( coord_y(j) < cg%coord_y(jj) - tolerance .AND. j < nyng ) 3235 j = j + 1 3117 3236 ENDDO 3118 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3119 jr = j 3120 DO WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy ) .AND.& 3121 ( j < nyng+1 ) ) 3122 j = j + 1 3123 jr = MIN( j, nyng ) 3124 ENDDO 3125 jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng ) 3126 jstart = jflv(jj) 3237 IF ( ABS( coord_y(j) - cg%coord_y(jj) ) < tolerance ) THEN 3238 j = jstart 3239 DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng ) 3240 j = j + 1 3241 ENDDO 3242 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3243 jfuv(jj) = jflv(jj) 3244 jstart = jflv(jj) 3245 ELSE 3246 j = jstart 3247 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) & 3248 .AND. ( j < nyng ) ) 3249 j = j + 1 3250 ENDDO 3251 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3252 jr = j 3253 DO WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy ) & 3254 .AND. ( j < nyng+1 ) ) 3255 j = j + 1 3256 jr = MIN( j, nyng ) 3257 ENDDO 3258 jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng ) 3259 jstart = jflv(jj) 3260 ENDIF 3261 !AH 3262 write(9,"('pmci_init_anterp_tophat, jj, jflv, jfuv: ', 3(i4,2x))") & 3263 jj, jflv(jj), jfuv(jj) 3264 flush(9) 3265 3127 3266 ENDDO 3128 3267 jflv(jcn) = nyng … … 3147 3286 jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng ) 3148 3287 jstart = jflo(jj) 3288 !AH 3289 write(9,"('pmci_init_anterp_tophat, ii, jflo, jfuo: ', 3(i4,2x))") & 3290 jj, jflo(jj), jfuo(jj) 3291 flush(9) 3292 3149 3293 ENDDO 3150 3294 jflo(jcn) = nyng … … 3152 3296 ! 3153 3297 !-- k-indices of w for each kk-index value 3298 !-- Note that anterpolation index limits are needed also for the top boundary 3299 !-- ghost cell level because of the reversibility correction in the interpolation. 3154 3300 kstart = 0 3155 3301 kflw(0) = 0 3156 3302 kfuw(0) = 0 3157 DO kk = 1, kctw 3303 tolerance = 0.000001_wp * dzw(1) 3304 DO kk = 1, cg%nz+1 3305 ! 3306 !-- In case the child and parent grid lines match in z 3307 !-- use only the local j,i-child-grid plane for the anterpolation, 3308 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3158 3309 k = kstart 3159 DO WHILE ( ( zw(k) < cg%zu(kk) ) .AND. ( k < nzt ))3310 DO WHILE ( zw(k) < cg%zw(kk) - tolerance .AND. k < nzt+1 ) 3160 3311 k = k + 1 3161 3312 ENDDO 3162 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3163 DO WHILE ( ( zw(k) <= cg%zu(kk+1) ) .AND. ( k < nzt+1 ) ) 3164 k = k + 1 3165 ENDDO 3166 kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 ) 3167 kstart = kflw(kk) 3313 IF ( ABS( zw(k) - cg%zw(kk) ) < tolerance ) THEN 3314 k = kstart 3315 DO WHILE ( ( zw(k) < cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 3316 k = k + 1 3317 ENDDO 3318 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3319 kfuw(kk) = kflw(kk) 3320 kstart = kflw(kk) 3321 ELSE 3322 k = kstart 3323 DO WHILE ( ( zw(k) < cg%zu(kk) ) .AND. ( k < nzt ) ) 3324 k = k + 1 3325 ENDDO 3326 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3327 IF ( kk+1 <= cg%nz+1 ) THEN 3328 DO WHILE ( ( zw(k) <= cg%zu(kk+1) ) .AND. ( k < nzt+1 ) ) 3329 k = k + 1 3330 IF ( k > nzt + 1 ) EXIT ! This EXIT is to prevent zu(k) from flowing over. 3331 ENDDO 3332 kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 ) 3333 ELSE 3334 kfuw(kk) = kflw(kk) 3335 ENDIF 3336 kstart = kflw(kk) 3337 ENDIF 3338 !AH 3339 write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 3(i4,2x))") & 3340 kk, kflw(kk), kfuw(kk) 3341 flush(9) 3342 3168 3343 ENDDO 3169 3344 ! … … 3172 3347 kflo(0) = 0 3173 3348 kfuo(0) = 0 3174 DO kk = 1, kctu 3349 ! 3350 !-- Note that anterpolation index limits are needed also for the top boundary 3351 !-- ghost cell level because of the reversibility correction in the interpolation. 3352 !AH DO kk = 1, kcto+1 3353 DO kk = 1, cg%nz+1 3175 3354 k = kstart 3176 DO WHILE ( ( zu(k) < cg%zw(kk-1) ) .AND. ( k < nzt ) ) 3355 !AH DO WHILE ( ( zu(k) < cg%zw(kk-1) ) .AND. ( k < nzt ) ) 3356 !-- Note that this is an IMPORTANT correction for the reversibility correction 3357 DO WHILE ( ( zu(k) < cg%zw(kk-1) ) .AND. ( k <= nzt ) ) 3177 3358 k = k + 1 3178 3359 ENDDO 3179 3360 kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3180 DO WHILE ( ( zu(k) <= cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 3361 !AH DO WHILE ( ( zu(k) <= cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 3362 !-- Note that this is an IMPORTANT correction for the reversibility correction 3363 DO WHILE ( ( zu(k) <= cg%zw(kk) ) .AND. ( k <= nzt+1 ) ) 3181 3364 k = k + 1 3365 IF ( k > nzt + 1 ) EXIT ! This EXIT is to prevent zu(k) from flowing over. 3182 3366 ENDDO 3183 3367 kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 ) 3184 3368 kstart = kflo(kk) 3369 !AH 3370 write(9,"('init kflo, kfuo: ', 4(i3,2x), e12.5)") kk, kflo(kk), kfuo(kk), nzt, cg%zw(kk) 3371 flush(9) 3185 3372 ENDDO 3186 3373 ! 3187 !-- Precomputation of number of fine-grid nodes inside coarse-grid ij-faces.3188 !-- Note that ii, jj, and kk are coarse-grid indices.3374 !-- Precomputation of number of fine-grid nodes inside parent-grid cells. 3375 !-- Note that ii, jj, and kk are parent-grid indices. 3189 3376 !-- This information is needed in anterpolation. 3190 3377 DO ii = icl, icr 3191 3378 DO jj = jcs, jcn 3192 DO kk = 0, kctu 3379 !AH DO kk = 0, kcto+1 3380 DO kk = 0, cg%nz+1 3193 3381 ! 3194 3382 !-- u-component … … 3223 3411 ENDDO 3224 3412 3225 DO kk = 0, kctw 3413 !AH DO kk = 0, kctw+1 3414 DO kk = 0, cg%nz+1 3226 3415 ! 3227 3416 !-- w-component … … 3259 3448 ENDIF 3260 3449 ENDDO 3261 3262 3450 3263 3451 DO jj = jcs, jcn … … 3275 3463 ENDIF 3276 3464 3277 ALLOCATE( fraz(0:kct u) )3278 DO kk = 0, kct u3465 ALLOCATE( fraz(0:kcto) ) 3466 DO kk = 0, kcto 3279 3467 zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4 3280 3468 fraz(kk) = zeta / ( 1.0_wp + zeta ) … … 3462 3650 coord_x(i) = lower_left_coord_x + i * dx 3463 3651 ENDDO 3464 3652 3465 3653 DO j = -nbgp, ny + nbgp 3466 3654 coord_y(j) = lower_left_coord_y + j * dy … … 3472 3660 3473 3661 3474 3475 3662 SUBROUTINE pmci_set_array_pointer( name, child_id, nz_cl, n ) 3476 3663 … … 3487 3674 3488 3675 REAL(wp), POINTER, DIMENSION(:,:) :: p_2d !< 3676 REAL(wp), POINTER, DIMENSION(:,:) :: p_2d_sec !< 3489 3677 REAL(wp), POINTER, DIMENSION(:,:,:) :: p_3d !< 3490 3678 REAL(wp), POINTER, DIMENSION(:,:,:) :: p_3d_sec !< … … 3584 3772 3585 3773 3774 3586 3775 INTEGER FUNCTION get_number_of_childs () 3587 3776 … … 3614 3803 3615 3804 END FUNCTION get_childid 3805 3616 3806 3617 3807 … … 3641 3831 END SUBROUTINE get_child_edges 3642 3832 3833 3834 3643 3835 SUBROUTINE get_child_gridspacing (m, dx,dy,dz) 3644 3836 … … 3656 3848 END SUBROUTINE get_child_gridspacing 3657 3849 3658 SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc,n ) 3850 3851 3852 SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc, n ) 3659 3853 3660 3854 IMPLICIT NONE … … 3666 3860 INTEGER(iwp), INTENT(IN) :: je !< 3667 3861 INTEGER(iwp), INTENT(IN) :: js !< 3668 INTEGER(iwp), INTENT(IN) :: nzc !< Note that nzc is cg%nz3862 INTEGER(iwp), INTENT(IN) :: nzc !< nzc is cg%nz, but note that cg%nz is not the original nz of parent, but the highest parent-grid level needed for nesting. 3669 3863 3670 3864 INTEGER(iwp), INTENT(IN), OPTIONAL :: n !< number of chemical species … … 3819 4013 !-- The interpolation. 3820 4014 CALL pmci_interp_tril_all ( u, uc, icu, jco, kco, r1xu, r2xu, r1yo, & 3821 r2yo, r1zo, r2zo, 'u' ) 4015 r2yo, r1zo, r2zo, kcto, iflu, ifuu, & 4016 jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' ) 3822 4017 CALL pmci_interp_tril_all ( v, vc, ico, jcv, kco, r1xo, r2xo, r1yv, & 3823 r2yv, r1zo, r2zo, 'v' ) 4018 r2yv, r1zo, r2zo, kcto, iflo, ifuo, & 4019 jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' ) 3824 4020 CALL pmci_interp_tril_all ( w, wc, ico, jco, kcw, r1xo, r2xo, r1yo, & 3825 r2yo, r1zw, r2zw, 'w' ) 4021 r2yo, r1zw, r2zw, kctw, iflo, ifuo, & 4022 jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' ) 3826 4023 3827 4024 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & … … 3829 4026 .NOT. constant_diffusion ) ) THEN 3830 4027 CALL pmci_interp_tril_all ( e, ec, ico, jco, kco, r1xo, r2xo, r1yo, & 3831 r2yo, r1zo, r2zo, 'e' ) 4028 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 4029 jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' ) 3832 4030 ENDIF 3833 4031 3834 4032 IF ( rans_mode_parent .AND. rans_mode .AND. rans_tke_e ) THEN 3835 4033 CALL pmci_interp_tril_all ( diss, dissc, ico, jco, kco, r1xo, r2xo,& 3836 r1yo, r2yo, r1zo, r2zo, 's' ) 4034 r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,& 4035 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 3837 4036 ENDIF 3838 4037 3839 4038 IF ( .NOT. neutral ) THEN 3840 4039 CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, & 3841 r1yo, r2yo, r1zo, r2zo, 's' ) 4040 r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,& 4041 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 3842 4042 ENDIF 3843 4043 … … 3845 4045 3846 4046 CALL pmci_interp_tril_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, & 3847 r2yo, r1zo, r2zo, 's' ) 4047 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 4048 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 3848 4049 3849 4050 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN 3850 4051 CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo, & 3851 r1yo, r2yo, r1zo, r2zo, 's' ) 4052 r1yo, r2yo, r1zo, r2zo, kcto, & 4053 iflo, ifuo, jflo, jfuo, kflo, kfuo, & 4054 ijkfc_s, 's' ) 3852 4055 CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo, & 3853 r1yo, r2yo, r1zo, r2zo, 's' ) 4056 r1yo, r2yo, r1zo, r2zo, kcto, & 4057 iflo, ifuo, jflo, jfuo, kflo, kfuo, & 4058 ijkfc_s, 's' ) 3854 4059 ENDIF 3855 4060 3856 4061 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN 3857 4062 CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo, & 3858 r1yo, r2yo, r1zo, r2zo, 's' ) 4063 r1yo, r2yo, r1zo, r2zo, kcto, & 4064 iflo, ifuo, jflo, jfuo, kflo, kfuo, & 4065 ijkfc_s, 's' ) 3859 4066 CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo, & 3860 r1yo, r2yo, r1zo, r2zo, 's' ) 4067 r1yo, r2yo, r1zo, r2zo, kcto, & 4068 iflo, ifuo, jflo, jfuo, kflo, kfuo, & 4069 ijkfc_s, 's' ) 3861 4070 ENDIF 3862 4071 … … 3865 4074 IF ( passive_scalar ) THEN 3866 4075 CALL pmci_interp_tril_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo, & 3867 r2yo, r1zo, r2zo, 's' ) 4076 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 4077 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 3868 4078 ENDIF 3869 4079 … … 3873 4083 chem_spec_c(:,:,:,n), & 3874 4084 ico, jco, kco, r1xo, r2xo, r1yo, & 3875 r2yo, r1zo, r2zo, 's' ) 4085 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 4086 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 3876 4087 ENDDO 3877 4088 ENDIF … … 3912 4123 3913 4124 SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, & 3914 r1z, r2z, var ) 4125 r1z, r2z, kct, ifl, ifu, jfl, jfu, & 4126 kfl, kfu, ijkfc, var ) 3915 4127 ! 3916 4128 !-- Interpolation of the internal values for the child-domain initialization … … 3925 4137 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !< 3926 4138 4139 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !< 4140 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc !< 4141 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !< 4142 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !< 4143 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !< 4144 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !< 4145 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !< 4146 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !< 4147 4148 INTEGER(iwp) :: kct 4149 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !< Indicates start index of child cells belonging to certain parent cell - x direction 4150 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !< Indicates end index of child cells belonging to certain parent cell - x direction 4151 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !< Indicates start index of child cells belonging to certain parent cell - y direction 4152 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !< Indicates start index of child cells belonging to certain parent cell - y direction 4153 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 4154 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 4155 INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 4156 3927 4157 INTEGER(iwp) :: i !< 3928 4158 INTEGER(iwp) :: ib !< 3929 4159 INTEGER(iwp) :: ie !< 4160 INTEGER(iwp) :: ijk !< 3930 4161 INTEGER(iwp) :: j !< 3931 4162 INTEGER(iwp) :: jb !< … … 3939 4170 INTEGER(iwp) :: m !< 3940 4171 INTEGER(iwp) :: n !< 3941 3942 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !< 3943 REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc !< 3944 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x !< 3945 REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x !< 3946 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y !< 3947 REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y !< 3948 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z !< 3949 REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z !< 3950 4172 INTEGER(iwp) :: var_flag !< 4173 4174 REAL(wp) :: cellsum !< 4175 REAL(wp) :: cellsumd !< 3951 4176 REAL(wp) :: fk !< 3952 4177 REAL(wp) :: fkj !< … … 3957 4182 REAL(wp) :: logratio !< 3958 4183 REAL(wp) :: logzuc1 !< 4184 REAL(wp) :: rcorr !< 4185 REAL(wp) :: rcorr_ijk !< 3959 4186 REAL(wp) :: zuc1 !< 3960 4187 REAL(wp) :: z0_topo !< roughness at vertical walls … … 3965 4192 jb = nys 3966 4193 je = nyn 4194 kb = 0 3967 4195 IF ( nesting_mode /= 'vertical' ) THEN 3968 4196 IF ( bc_dirichlet_l ) THEN … … 3989 4217 ENDIF 3990 4218 ENDIF 4219 4220 IF ( var == 'u' ) THEN 4221 var_flag = 1 4222 ELSEIF ( var == 'v' ) THEN 4223 var_flag = 2 4224 ELSEIF ( var == 'w' ) THEN 4225 var_flag = 3 4226 ELSE 4227 var_flag = 0 4228 ENDIF 3991 4229 ! 3992 4230 !-- Trilinear interpolation. … … 3997 4235 !-- topography top at grid point (j,i) in order to not overwrite 3998 4236 !-- the bottom BC. 3999 kb = get_topography_top_index_ji( j, i, TRIM ( var ) ) + 1 4237 ! kb = get_topography_top_index_ji( j, i, TRIM ( var ) ) + 1 4000 4238 DO k = kb, nzt + 1 4001 4239 l = ic(i) … … 4062 4300 4063 4301 ENDIF 4302 ! 4303 !-- Apply the reversibility correction. 4304 DO l = icl, icr 4305 DO m = jcs, jcn 4306 DO n = 0, kct+1 4307 ijk = 1 4308 cellsum = 0.0_wp 4309 cellsumd = 0.0_wp 4310 ! 4311 !-- Note that the index name i must not be used here as a loop 4312 !-- index name since i is the constant boundary index, hence 4313 !-- the name ia. 4314 DO i = ifl(l), ifu(l) 4315 DO j = jfl(m), jfu(m) 4316 DO k = kfl(n), kfu(n) 4317 cellsum = cellsum + MERGE( f(k,j,i), 0.0_wp, & 4318 BTEST( wall_flags_0(k,j,i), var_flag ) ) 4319 celltmpd(ijk) = ABS( fc(n,m,l) - f(k,j,i) ) 4320 cellsumd = cellsumd + MERGE( celltmpd(ijk), & 4321 0.0_wp, BTEST( wall_flags_0(k,j,i), var_flag ) ) 4322 ijk = ijk + 1 4323 ENDDO 4324 ENDDO 4325 ENDDO 4326 4327 IF ( ijkfc(n,m,l) /= 0 ) THEN 4328 cellsum = cellsum / REAL( ijkfc(n,m,l), KIND=wp ) 4329 rcorr = fc(n,m,l) - cellsum 4330 cellsumd = cellsumd / REAL( ijkfc(n,m,l), KIND=wp ) 4331 ELSE 4332 cellsum = 0.0_wp 4333 rcorr = 0.0_wp 4334 cellsumd = 1.0_wp 4335 celltmpd = 1.0_wp 4336 ENDIF 4337 ! 4338 !-- Distribute the correction term to the child nodes according to 4339 !-- their relative difference to the parent value such that the 4340 !-- node with the largest difference gets the largest share of the 4341 !-- correction. The distribution is skipped if rcorr is negligibly 4342 !-- small in order to avoid division by zero. 4343 IF ( ABS(rcorr) < 0.000001_wp ) THEN 4344 cellsumd = 1.0_wp 4345 celltmpd = 1.0_wp 4346 ENDIF 4347 4348 ijk = 1 4349 DO i = ifl(l), ifu(l) 4350 DO j = jfl(m), jfu(m) 4351 DO k = kfl(n), kfu(n) 4352 rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 4353 f(k,j,i) = f(k,j,i) + rcorr_ijk 4354 ijk = ijk + 1 4355 ENDDO 4356 ENDDO 4357 ENDDO 4358 4359 ENDDO ! n 4360 ENDDO ! m 4361 ENDDO ! l 4064 4362 4065 4363 END SUBROUTINE pmci_interp_tril_all … … 4165 4463 4166 4464 4167 SUBROUTINE pmci_ensure_nest_mass_conservation4168 4169 !4170 !-- Adjust the volume-flow rate through the top boundary so that the net volume4171 !-- flow through all boundaries of the current nest domain becomes zero.4172 IMPLICIT NONE4173 4174 INTEGER(iwp) :: i !<4175 INTEGER(iwp) :: ierr !<4176 INTEGER(iwp) :: j !<4177 INTEGER(iwp) :: k !<4178 4179 REAL(wp) :: dxdy !<4180 REAL(wp) :: innor !<4181 REAL(wp) :: w_lt !<4182 REAL(wp), DIMENSION(1:3) :: volume_flow_l !<4183 4184 !4185 !-- Sum up the volume flow through the left/right boundaries4186 volume_flow(1) = 0.0_wp4187 volume_flow_l(1) = 0.0_wp4188 4189 IF ( bc_dirichlet_l ) THEN4190 i = 04191 innor = dy4192 DO j = nys, nyn4193 DO k = nzb+1, nzt4194 volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k) &4195 * MERGE( 1.0_wp, 0.0_wp, &4196 BTEST( wall_flags_0(k,j,i), 1 ) )4197 ENDDO4198 ENDDO4199 ENDIF4200 4201 IF ( bc_dirichlet_r ) THEN4202 i = nx + 14203 innor = -dy4204 DO j = nys, nyn4205 DO k = nzb+1, nzt4206 volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k) &4207 * MERGE( 1.0_wp, 0.0_wp, &4208 BTEST( wall_flags_0(k,j,i), 1 ) )4209 ENDDO4210 ENDDO4211 ENDIF4212 4213 #if defined( __parallel )4214 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )4215 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &4216 MPI_SUM, comm2d, ierr )4217 #else4218 volume_flow(1) = volume_flow_l(1)4219 #endif4220 4221 !4222 !-- Sum up the volume flow through the south/north boundaries4223 volume_flow(2) = 0.0_wp4224 volume_flow_l(2) = 0.0_wp4225 4226 IF ( bc_dirichlet_s ) THEN4227 j = 04228 innor = dx4229 DO i = nxl, nxr4230 DO k = nzb+1, nzt4231 volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k) &4232 * MERGE( 1.0_wp, 0.0_wp, &4233 BTEST( wall_flags_0(k,j,i), 2 ) )4234 ENDDO4235 ENDDO4236 ENDIF4237 4238 IF ( bc_dirichlet_n ) THEN4239 j = ny + 14240 innor = -dx4241 DO i = nxl, nxr4242 DO k = nzb+1, nzt4243 volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k) &4244 * MERGE( 1.0_wp, 0.0_wp, &4245 BTEST( wall_flags_0(k,j,i), 2 ) )4246 ENDDO4247 ENDDO4248 ENDIF4249 4250 #if defined( __parallel )4251 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )4252 CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &4253 MPI_SUM, comm2d, ierr )4254 #else4255 volume_flow(2) = volume_flow_l(2)4256 #endif4257 4258 !4259 !-- Sum up the volume flow through the top boundary4260 volume_flow(3) = 0.0_wp4261 volume_flow_l(3) = 0.0_wp4262 dxdy = dx * dy4263 k = nzt4264 DO i = nxl, nxr4265 DO j = nys, nyn4266 volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy4267 ENDDO4268 ENDDO4269 4270 #if defined( __parallel )4271 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )4272 CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL, &4273 MPI_SUM, comm2d, ierr )4274 #else4275 volume_flow(3) = volume_flow_l(3)4276 #endif4277 4278 !4279 !-- Correct the top-boundary value of w4280 w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t4281 DO i = nxl, nxr4282 DO j = nys, nyn4283 DO k = nzt, nzt + 14284 w(k,j,i) = w(k,j,i) + w_lt4285 ENDDO4286 ENDDO4287 ENDDO4288 4289 END SUBROUTINE pmci_ensure_nest_mass_conservation4290 4291 4292 4293 4465 SUBROUTINE pmci_synchronize 4294 4466 … … 4409 4581 DO m = 1, SIZE( pmc_parent_for_child ) - 1 4410 4582 child_id = pmc_parent_for_child(m) 4411 4412 4583 IF ( direction == parent_to_child ) THEN 4413 4584 CALL cpu_log( log_point_s(71), 'pmc parent send', 'start' ) … … 4484 4655 4485 4656 IF ( direction == parent_to_child ) THEN 4486 4657 4487 4658 CALL cpu_log( log_point_s(73), 'pmc child recv', 'start' ) 4488 4659 CALL pmc_c_getbuffer( ) … … 4492 4663 CALL pmci_interpolation 4493 4664 CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' ) 4494 4665 4495 4666 ELSE 4496 4667 ! … … 4517 4688 IMPLICIT NONE 4518 4689 4690 INTEGER(iwp) :: ibgp !< index running over the nbgp boundary ghost points in i-direction 4691 INTEGER(iwp) :: jbgp !< index running over the nbgp boundary ghost points in j-direction 4519 4692 INTEGER(iwp) :: n !< running index for number of chemical species 4520 4693 … … 4527 4700 !-- Left border pe: 4528 4701 IF ( bc_dirichlet_l ) THEN 4529 4530 CALL pmci_interp_tril_lr( u, uc, icu, jco, kco, r1xu, r2xu, & 4531 r1yo, r2yo, r1zo, r2zo, & 4532 logc_u_l, logc_ratio_u_l, & 4533 logc_kbounds_u_l, & 4534 nzt_topo_nestbc_l, 'l', 'u' ) 4535 4536 CALL pmci_interp_tril_lr( v, vc, ico, jcv, kco, r1xo, r2xo, & 4537 r1yv, r2yv, r1zo, r2zo, & 4538 logc_v_l, logc_ratio_v_l, & 4539 logc_kbounds_v_l, & 4540 nzt_topo_nestbc_l, 'l', 'v' ) 4541 4542 CALL pmci_interp_tril_lr( w, wc, ico, jco, kcw, r1xo, r2xo, & 4543 r1yo, r2yo, r1zw, r2zw, & 4544 logc_w_l, logc_ratio_w_l, & 4545 logc_kbounds_w_l, & 4546 nzt_topo_nestbc_l, 'l', 'w' ) 4702 4703 CALL pmci_interp_tril_lr( u, uc, icu, jco, kco, r1xu, r2xu, & 4704 r1yo, r2yo, r1zo, r2zo, & 4705 logc_u_l, logc_ratio_u_l, & 4706 logc_kbounds_u_l, nzt_topo_nestbc_l, & 4707 kcto, iflu, ifuu, jflo, jfuo, kflo, & 4708 kfuo, ijkfc_u, 'l', 'u' ) 4709 4710 CALL pmci_interp_tril_lr( v, vc, ico, jcv, kco, r1xo, r2xo, & 4711 r1yv, r2yv, r1zo, r2zo, & 4712 logc_v_l, logc_ratio_v_l, & 4713 logc_kbounds_v_l, nzt_topo_nestbc_l, & 4714 kcto, iflo, ifuo, jflv, jfuv, kflo, & 4715 kfuo, ijkfc_v, 'l', 'v' ) 4716 4717 CALL pmci_interp_tril_lr( w, wc, ico, jco, kcw, r1xo, r2xo, & 4718 r1yo, r2yo, r1zw, r2zw, & 4719 logc_w_l, logc_ratio_w_l, & 4720 logc_kbounds_w_l, nzt_topo_nestbc_l, & 4721 kctw, iflo, ifuo, jflo, jfuo, kflw, & 4722 kfuw, ijkfc_w, 'l', 'w' ) 4547 4723 4548 4724 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & 4549 4725 ( .NOT. rans_mode_parent .AND. .NOT. rans_mode .AND. & 4550 4726 .NOT. constant_diffusion ) ) THEN 4551 CALL pmci_interp_tril_lr( e, ec, ico, jco, kco, r1xo, r2xo, & 4552 r1yo, r2yo, r1zo, r2zo, & 4553 logc_w_l, logc_ratio_w_l, & 4554 logc_kbounds_w_l, & 4555 nzt_topo_nestbc_l, 'l', 'e' ) 4727 ! CALL pmci_interp_tril_lr( e, ec, ico, jco, kco, r1xo, r2xo, & 4728 ! r1yo, r2yo, r1zo, r2zo, & 4729 ! logc_w_l, logc_ratio_w_l, & 4730 ! logc_kbounds_w_l, nzt_topo_nestbc_l, & 4731 ! kcto, iflo, ifuo, jflo, jfuo, kflo, & 4732 ! kfuo, ijkfc_s, 'l', 'e' ) 4733 ! 4734 !-- Interpolation of e is replaced by the Neumann condition. 4735 DO ibgp = -nbgp, -1 4736 e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,0) 4737 ENDDO 4738 4556 4739 ENDIF 4557 4740 … … 4560 4743 r2xo, r1yo, r2yo, r1zo, r2zo, & 4561 4744 logc_w_l, logc_ratio_w_l, & 4562 logc_kbounds_w_l, & 4563 nzt_topo_nestbc_l, 'l', 's' ) 4745 logc_kbounds_w_l, nzt_topo_nestbc_l, & 4746 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4747 kfuo, ijkfc_s, 'l', 's' ) 4564 4748 ENDIF 4565 4749 … … 4568 4752 r1yo, r2yo, r1zo, r2zo, & 4569 4753 logc_w_l, logc_ratio_w_l, & 4570 logc_kbounds_w_l, & 4571 nzt_topo_nestbc_l, 'l', 's' ) 4754 logc_kbounds_w_l, nzt_topo_nestbc_l, & 4755 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4756 kfuo, ijkfc_s, 'l', 's' ) 4572 4757 ENDIF 4573 4758 … … 4577 4762 r1yo, r2yo, r1zo, r2zo, & 4578 4763 logc_w_l, logc_ratio_w_l, & 4579 logc_kbounds_w_l, & 4580 nzt_topo_nestbc_l, 'l', 's' ) 4764 logc_kbounds_w_l, nzt_topo_nestbc_l, & 4765 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4766 kfuo, ijkfc_s, 'l', 's' ) 4581 4767 4582 4768 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN … … 4585 4771 logc_w_l, logc_ratio_w_l, & 4586 4772 logc_kbounds_w_l, & 4587 nzt_topo_nestbc_l, 'l', 's' ) 4773 nzt_topo_nestbc_l, & 4774 kcto, iflo, ifuo, jflo, jfuo, & 4775 kflo, kfuo, ijkfc_s, 'l', 's' ) 4588 4776 4589 4777 CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, & … … 4591 4779 logc_w_l, logc_ratio_w_l, & 4592 4780 logc_kbounds_w_l, & 4593 nzt_topo_nestbc_l, 'l', 's' ) 4781 nzt_topo_nestbc_l, & 4782 kcto, iflo, ifuo, jflo, jfuo, & 4783 kflo, kfuo, ijkfc_s, 'l', 's' ) 4594 4784 ENDIF 4595 4785 … … 4599 4789 logc_w_l, logc_ratio_w_l, & 4600 4790 logc_kbounds_w_l, & 4601 nzt_topo_nestbc_l, 'l', 's' ) 4791 nzt_topo_nestbc_l, & 4792 kcto, iflo, ifuo, jflo, jfuo, & 4793 kflo, kfuo, ijkfc_s, 'l', 's' ) 4602 4794 4603 4795 CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo, & 4604 4796 r2xo, r1yo, r2yo, r1zo, r2zo, & 4605 4797 logc_w_l, logc_ratio_w_l, & 4606 logc_kbounds_w_l, & 4607 nzt_topo_nestbc_l, 'l', 's' ) 4798 logc_kbounds_w_l, & 4799 nzt_topo_nestbc_l, & 4800 kcto, iflo, ifuo, jflo, jfuo, & 4801 kflo, kfuo, ijkfc_s, 'l', 's' ) 4608 4802 ENDIF 4609 4803 … … 4615 4809 logc_w_l, logc_ratio_w_l, & 4616 4810 logc_kbounds_w_l, & 4617 nzt_topo_nestbc_l, 'l', 's' ) 4811 nzt_topo_nestbc_l, & 4812 kcto, iflo, ifuo, jflo, jfuo, & 4813 kflo, kfuo, ijkfc_s, 'l', 's' ) 4618 4814 ENDIF 4619 4815 … … 4626 4822 logc_w_l, logc_ratio_w_l, & 4627 4823 logc_kbounds_w_l, & 4628 nzt_topo_nestbc_l, 'l', 's' ) 4824 nzt_topo_nestbc_l, & 4825 kcto, iflo, ifuo, jflo, jfuo, & 4826 kflo, kfuo, ijkfc_s, 'l', 's' ) 4629 4827 ENDDO 4630 4828 ENDIF … … 4638 4836 r1yo, r2yo, r1zo, r2zo, & 4639 4837 logc_u_r, logc_ratio_u_r, & 4640 logc_kbounds_u_r, & 4641 nzt_topo_nestbc_r, 'r', 'u' ) 4838 logc_kbounds_u_r, nzt_topo_nestbc_r, & 4839 kcto, iflu, ifuu, jflo, jfuo, kflo, & 4840 kfuo, ijkfc_u, 'r', 'u' ) 4642 4841 4643 4842 CALL pmci_interp_tril_lr( v, vc, ico, jcv, kco, r1xo, r2xo, & 4644 4843 r1yv, r2yv, r1zo, r2zo, & 4645 4844 logc_v_r, logc_ratio_v_r, & 4646 logc_kbounds_v_r, & 4647 nzt_topo_nestbc_r, 'r', 'v' ) 4845 logc_kbounds_v_r, nzt_topo_nestbc_r, & 4846 kcto, iflo, ifuo, jflv, jfuv, kflo, & 4847 kfuo, ijkfc_v, 'r', 'v' ) 4648 4848 4649 4849 CALL pmci_interp_tril_lr( w, wc, ico, jco, kcw, r1xo, r2xo, & 4650 4850 r1yo, r2yo, r1zw, r2zw, & 4651 4851 logc_w_r, logc_ratio_w_r, & 4652 logc_kbounds_w_r, & 4653 nzt_topo_nestbc_r, 'r', 'w' ) 4852 logc_kbounds_w_r, nzt_topo_nestbc_r, & 4853 kctw, iflo, ifuo, jflo, jfuo, kflw, & 4854 kfuw, ijkfc_w, 'r', 'w' ) 4654 4855 4655 4856 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & 4656 4857 ( .NOT. rans_mode_parent .AND. .NOT. rans_mode .AND. & 4657 4858 .NOT. constant_diffusion ) ) THEN 4658 CALL pmci_interp_tril_lr( e, ec, ico, jco, kco, r1xo, r2xo, & 4659 r1yo,r2yo, r1zo, r2zo, & 4660 logc_w_r, logc_ratio_w_r, & 4661 logc_kbounds_w_r, & 4662 nzt_topo_nestbc_r, 'r', 'e' ) 4663 4859 ! CALL pmci_interp_tril_lr( e, ec, ico, jco, kco, r1xo, r2xo, & 4860 ! r1yo,r2yo, r1zo, r2zo, & 4861 ! logc_w_r, logc_ratio_w_r, & 4862 ! logc_kbounds_w_r, nzt_topo_nestbc_r, & 4863 ! kcto, iflo, ifuo, jflo, jfuo, kflo, & 4864 ! kfuo, ijkfc_s, 'r', 'e' ) 4865 ! 4866 !-- Interpolation of e is replaced by the Neumann condition. 4867 DO ibgp = nx+1, nx+nbgp 4868 e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,nx) 4869 ENDDO 4664 4870 ENDIF 4665 4871 … … 4668 4874 r2xo, r1yo,r2yo, r1zo, r2zo, & 4669 4875 logc_w_r, logc_ratio_w_r, & 4670 logc_kbounds_w_r, & 4671 nzt_topo_nestbc_r, 'r', 's' ) 4876 logc_kbounds_w_r, nzt_topo_nestbc_r, & 4877 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4878 kfuo, ijkfc_s, 'r', 's' ) 4672 4879 4673 4880 ENDIF … … 4677 4884 r1yo, r2yo, r1zo, r2zo, & 4678 4885 logc_w_r, logc_ratio_w_r, & 4679 logc_kbounds_w_r, & 4680 nzt_topo_nestbc_r, 'r', 's' ) 4886 logc_kbounds_w_r, nzt_topo_nestbc_r, & 4887 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4888 kfuo, ijkfc_s, 'r', 's' ) 4681 4889 ENDIF 4682 4890 … … 4685 4893 r1yo, r2yo, r1zo, r2zo, & 4686 4894 logc_w_r, logc_ratio_w_r, & 4687 logc_kbounds_w_r, & 4688 nzt_topo_nestbc_r, 'r', 's' ) 4895 logc_kbounds_w_r, nzt_topo_nestbc_r, & 4896 kcto, iflo, ifuo, jflo, jfuo, kflo, & 4897 kfuo, ijkfc_s, 'r', 's' ) 4689 4898 4690 4899 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN … … 4694 4903 logc_w_r, logc_ratio_w_r, & 4695 4904 logc_kbounds_w_r, & 4696 nzt_topo_nestbc_r, 'r', 's' ) 4905 nzt_topo_nestbc_r, & 4906 kcto, iflo, ifuo, jflo, jfuo, & 4907 kflo, kfuo, ijkfc_s, 'r', 's' ) 4697 4908 4698 4909 CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, & … … 4700 4911 logc_w_r, logc_ratio_w_r, & 4701 4912 logc_kbounds_w_r, & 4702 nzt_topo_nestbc_r, 'r', 's' ) 4703 4913 nzt_topo_nestbc_r, & 4914 kcto, iflo, ifuo, jflo, jfuo, & 4915 kflo, kfuo, ijkfc_s, 'r', 's' ) 4704 4916 4705 4917 ENDIF … … 4712 4924 logc_w_r, logc_ratio_w_r, & 4713 4925 logc_kbounds_w_r, & 4714 nzt_topo_nestbc_r, 'r', 's' ) 4926 nzt_topo_nestbc_r, & 4927 kcto, iflo, ifuo, jflo, jfuo, & 4928 kflo, kfuo, ijkfc_s, & 4929 'r', 's' ) 4715 4930 4716 4931 CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo, & … … 4718 4933 logc_w_r, logc_ratio_w_r, & 4719 4934 logc_kbounds_w_r, & 4720 nzt_topo_nestbc_r, 'r', 's' ) 4935 nzt_topo_nestbc_r, & 4936 kcto, iflo, ifuo, jflo, jfuo, & 4937 kflo, kfuo, ijkfc_s, 'r', 's' ) 4721 4938 4722 4939 ENDIF … … 4729 4946 logc_w_r, logc_ratio_w_r, & 4730 4947 logc_kbounds_w_r, & 4731 nzt_topo_nestbc_r, 'r', 's' ) 4948 nzt_topo_nestbc_r, & 4949 kcto, iflo, ifuo, jflo, jfuo, & 4950 kflo, kfuo, ijkfc_s, 'r', 's' ) 4732 4951 ENDIF 4733 4952 … … 4740 4959 logc_w_r, logc_ratio_w_r, & 4741 4960 logc_kbounds_w_r, & 4742 nzt_topo_nestbc_r, 'r', 's' ) 4961 nzt_topo_nestbc_r, & 4962 kcto, iflo, ifuo, jflo, jfuo, & 4963 kflo, kfuo, ijkfc_s, 'r', 's' ) 4743 4964 ENDDO 4744 4965 ENDIF … … 4751 4972 r1yo, r2yo, r1zo, r2zo, & 4752 4973 logc_u_s, logc_ratio_u_s, & 4753 logc_kbounds_u_s, & 4754 nzt_topo_nestbc_s, 's', 'u' ) 4974 logc_kbounds_u_s, nzt_topo_nestbc_s, & 4975 kcto, iflu, ifuu, jflo, jfuo, kflo, & 4976 kfuo, ijkfc_u, 's', 'u' ) 4755 4977 4756 4978 CALL pmci_interp_tril_sn( v, vc, ico, jcv, kco, r1xo, r2xo, & 4757 4979 r1yv, r2yv, r1zo, r2zo, & 4758 4980 logc_v_s, logc_ratio_v_s, & 4759 logc_kbounds_v_s, & 4760 nzt_topo_nestbc_s, 's', 'v' ) 4981 logc_kbounds_v_s, nzt_topo_nestbc_s, & 4982 kcto, iflo, ifuo, jflv, jfuv, kflo, & 4983 kfuo, ijkfc_v, 's', 'v' ) 4761 4984 4762 4985 CALL pmci_interp_tril_sn( w, wc, ico, jco, kcw, r1xo, r2xo, & 4763 4986 r1yo, r2yo, r1zw, r2zw, & 4764 4987 logc_w_s, logc_ratio_w_s, & 4765 logc_kbounds_w_s, & 4766 nzt_topo_nestbc_s, 's','w' ) 4988 logc_kbounds_w_s, nzt_topo_nestbc_s, & 4989 kctw, iflo, ifuo, jflo, jfuo, kflw, & 4990 kfuw, ijkfc_w, 's','w' ) 4767 4991 4768 4992 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & 4769 4993 ( .NOT. rans_mode_parent .AND. .NOT. rans_mode .AND. & 4770 4994 .NOT. constant_diffusion ) ) THEN 4771 CALL pmci_interp_tril_sn( e, ec, ico, jco, kco, r1xo, r2xo, & 4772 r1yo, r2yo, r1zo, r2zo, & 4773 logc_w_s, logc_ratio_w_s, & 4774 logc_kbounds_w_s, & 4775 nzt_topo_nestbc_s, 's', 'e' ) 4776 4995 ! CALL pmci_interp_tril_sn( e, ec, ico, jco, kco, r1xo, r2xo, & 4996 ! r1yo, r2yo, r1zo, r2zo, & 4997 ! logc_w_s, logc_ratio_w_s, & 4998 ! logc_kbounds_w_s, nzt_topo_nestbc_s, & 4999 ! kcto, iflo, ifuo, jflo, jfuo, kflo, & 5000 ! kfuo, ijkfc_s, 's', 'e' ) 5001 ! 5002 !-- Interpolation of e is replaced by the Neumann condition. 5003 DO jbgp = -nbgp, -1 5004 e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,0,nxl:nxr) 5005 ENDDO 4777 5006 ENDIF 4778 5007 … … 4781 5010 r2xo, r1yo, r2yo, r1zo, r2zo, & 4782 5011 logc_w_s, logc_ratio_w_s, & 4783 logc_kbounds_w_s, 4784 nzt_topo_nestbc_s, 's', 's' )4785 5012 logc_kbounds_w_s, nzt_topo_nestbc_s, & 5013 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5014 kfuo, ijkfc_s, 's', 's' ) 4786 5015 ENDIF 4787 5016 … … 4790 5019 r1yo, r2yo, r1zo, r2zo, & 4791 5020 logc_w_s, logc_ratio_w_s, & 4792 logc_kbounds_w_s, & 4793 nzt_topo_nestbc_s, 's', 's' ) 5021 logc_kbounds_w_s, nzt_topo_nestbc_s, & 5022 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5023 kfuo, ijkfc_s, 's', 's' ) 4794 5024 ENDIF 4795 5025 … … 4798 5028 r1yo,r2yo, r1zo, r2zo, & 4799 5029 logc_w_s, logc_ratio_w_s, & 4800 logc_kbounds_w_s, & 4801 nzt_topo_nestbc_s, 's', 's' ) 5030 logc_kbounds_w_s, nzt_topo_nestbc_s, & 5031 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5032 kfuo, ijkfc_s, 's', 's' ) 4802 5033 4803 5034 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN … … 4807 5038 logc_w_s, logc_ratio_w_s, & 4808 5039 logc_kbounds_w_s, & 4809 nzt_topo_nestbc_s, 's', 's' ) 5040 nzt_topo_nestbc_s, & 5041 kcto, iflo, ifuo, jflo, jfuo, & 5042 kflo, kfuo, ijkfc_s, 's', 's' ) 4810 5043 4811 5044 CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo, & … … 4813 5046 logc_w_s, logc_ratio_w_s, & 4814 5047 logc_kbounds_w_s, & 4815 nzt_topo_nestbc_s, 's', 's' ) 5048 nzt_topo_nestbc_s, & 5049 kcto, iflo, ifuo, jflo, jfuo, & 5050 kflo, kfuo, ijkfc_s, 's', 's' ) 4816 5051 4817 5052 ENDIF … … 4823 5058 logc_w_s, logc_ratio_w_s, & 4824 5059 logc_kbounds_w_s, & 4825 nzt_topo_nestbc_s, 's', 's' ) 5060 nzt_topo_nestbc_s, & 5061 kcto, iflo, ifuo, jflo, jfuo, & 5062 kflo, kfuo, ijkfc_s, 's', 's' ) 4826 5063 4827 5064 CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo, & … … 4829 5066 logc_w_s, logc_ratio_w_s, & 4830 5067 logc_kbounds_w_s, & 4831 nzt_topo_nestbc_s, 's', 's' ) 5068 nzt_topo_nestbc_s, & 5069 kcto, iflo, ifuo, jflo, jfuo, & 5070 kflo, kfuo, ijkfc_s, 's', 's' ) 4832 5071 4833 5072 ENDIF … … 4840 5079 logc_w_s, logc_ratio_w_s, & 4841 5080 logc_kbounds_w_s, & 4842 nzt_topo_nestbc_s, 's', 's' ) 5081 nzt_topo_nestbc_s, & 5082 kcto, iflo, ifuo, jflo, jfuo, & 5083 kflo, kfuo, ijkfc_s, 's', 's' ) 4843 5084 ENDIF 4844 5085 … … 4851 5092 logc_w_s, logc_ratio_w_s, & 4852 5093 logc_kbounds_w_s, & 4853 nzt_topo_nestbc_s, 's', 's' ) 5094 nzt_topo_nestbc_s, & 5095 kcto, iflo, ifuo, jflo, jfuo, & 5096 kflo, kfuo, ijkfc_s, 's', 's' ) 4854 5097 ENDDO 4855 5098 ENDIF … … 4862 5105 r1yo, r2yo, r1zo, r2zo, & 4863 5106 logc_u_n, logc_ratio_u_n, & 4864 logc_kbounds_u_n, & 4865 nzt_topo_nestbc_n, 'n', 'u' ) 5107 logc_kbounds_u_n, nzt_topo_nestbc_n, & 5108 kcto, iflu, ifuu, jflo, jfuo, kflo, & 5109 kfuo, ijkfc_u, 'n', 'u' ) 4866 5110 4867 5111 CALL pmci_interp_tril_sn( v, vc, ico, jcv, kco, r1xo, r2xo, & 4868 5112 r1yv, r2yv, r1zo, r2zo, & 4869 5113 logc_v_n, logc_ratio_v_n, & 4870 logc_kbounds_v_n, & 4871 nzt_topo_nestbc_n, 'n', 'v' ) 5114 logc_kbounds_v_n, nzt_topo_nestbc_n, & 5115 kcto, iflo, ifuo, jflv, jfuv, kflo, & 5116 kfuo, ijkfc_v, 'n', 'v' ) 4872 5117 4873 5118 CALL pmci_interp_tril_sn( w, wc, ico, jco, kcw, r1xo, r2xo, & 4874 5119 r1yo, r2yo, r1zw, r2zw, & 4875 5120 logc_w_n, logc_ratio_w_n, & 4876 logc_kbounds_w_n, & 4877 nzt_topo_nestbc_n, 'n', 'w' ) 5121 logc_kbounds_w_n, nzt_topo_nestbc_n, & 5122 kctw, iflo, ifuo, jflo, jfuo, kflw, & 5123 kfuw, ijkfc_w, 'n', 'w' ) 4878 5124 4879 5125 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & 4880 5126 ( .NOT. rans_mode_parent .AND. .NOT. rans_mode .AND. & 4881 5127 .NOT. constant_diffusion ) ) THEN 4882 CALL pmci_interp_tril_sn( e, ec, ico, jco, kco, r1xo, r2xo, & 4883 r1yo, r2yo, r1zo, r2zo, & 4884 logc_w_n, logc_ratio_w_n, & 4885 logc_kbounds_w_n, & 4886 nzt_topo_nestbc_n, 'n', 'e' ) 4887 5128 ! CALL pmci_interp_tril_sn( e, ec, ico, jco, kco, r1xo, r2xo, & 5129 ! r1yo, r2yo, r1zo, r2zo, & 5130 ! logc_w_n, logc_ratio_w_n, & 5131 ! logc_kbounds_w_n, nzt_topo_nestbc_n, & 5132 ! kcto, iflo, ifuo, jflo, jfuo, kflo, & 5133 ! kfuo, ijkfc_s, 'n', 'e' ) 5134 ! 5135 !-- Interpolation of e is replaced by the Neumann condition. 5136 DO jbgp = ny+1, ny+nbgp 5137 e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,ny,nxl:nxr) 5138 ENDDO 4888 5139 ENDIF 4889 5140 … … 4892 5143 r2xo, r1yo, r2yo, r1zo, r2zo, & 4893 5144 logc_w_n, logc_ratio_w_n, & 4894 logc_kbounds_w_n, & 4895 nzt_topo_nestbc_n, 'n', 's' ) 5145 logc_kbounds_w_n, nzt_topo_nestbc_n, & 5146 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5147 kfuo, ijkfc_s, 'n', 's' ) 4896 5148 4897 5149 ENDIF … … 4901 5153 r1yo, r2yo, r1zo, r2zo, & 4902 5154 logc_w_n, logc_ratio_w_n, & 4903 logc_kbounds_w_n, & 4904 nzt_topo_nestbc_n, 'n', 's' ) 5155 logc_kbounds_w_n, nzt_topo_nestbc_n, & 5156 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5157 kfuo, ijkfc_s, 'n', 's' ) 4905 5158 ENDIF 4906 5159 … … 4909 5162 r1yo, r2yo, r1zo, r2zo, & 4910 5163 logc_w_n, logc_ratio_w_n, & 4911 logc_kbounds_w_n, & 4912 nzt_topo_nestbc_n, 'n', 's' ) 5164 logc_kbounds_w_n, nzt_topo_nestbc_n, & 5165 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5166 kfuo, ijkfc_s, 'n', 's' ) 4913 5167 4914 5168 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN … … 4918 5172 logc_w_n, logc_ratio_w_n, & 4919 5173 logc_kbounds_w_n, & 4920 nzt_topo_nestbc_n, 'n', 's' ) 5174 nzt_topo_nestbc_n, & 5175 kcto, iflo, ifuo, jflo, jfuo, & 5176 kflo, kfuo, ijkfc_s, 'n', 's' ) 4921 5177 4922 5178 CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo, & … … 4924 5180 logc_u_n, logc_ratio_u_n, & 4925 5181 logc_kbounds_w_n, & 4926 nzt_topo_nestbc_n, 'n', 's' ) 5182 nzt_topo_nestbc_n, & 5183 kcto, iflo, ifuo, jflo, jfuo, & 5184 kflo, kfuo, ijkfc_s, 'n', 's' ) 4927 5185 4928 5186 ENDIF … … 4934 5192 logc_w_n, logc_ratio_w_n, & 4935 5193 logc_kbounds_w_n, & 4936 nzt_topo_nestbc_n, 'n', 's' ) 5194 nzt_topo_nestbc_n, & 5195 kcto, iflo, ifuo, jflo, jfuo, & 5196 kflo, kfuo, ijkfc_s, 'n', 's' ) 4937 5197 4938 5198 CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo, & … … 4940 5200 logc_w_n, logc_ratio_w_n, & 4941 5201 logc_kbounds_w_n, & 4942 nzt_topo_nestbc_n, 'n', 's' ) 5202 nzt_topo_nestbc_n, & 5203 kcto, iflo, ifuo, jflo, jfuo, & 5204 kflo, kfuo, ijkfc_s, 'n', 's' ) 4943 5205 4944 5206 ENDIF … … 4951 5213 logc_w_n, logc_ratio_w_n, & 4952 5214 logc_kbounds_w_n, & 4953 nzt_topo_nestbc_n, 'n', 's' ) 5215 nzt_topo_nestbc_n, & 5216 kcto, iflo, ifuo, jflo, jfuo, & 5217 kflo, kfuo, ijkfc_s, 'n', 's' ) 4954 5218 ENDIF 4955 5219 … … 4962 5226 logc_w_n, logc_ratio_w_n, & 4963 5227 logc_kbounds_w_n, & 4964 nzt_topo_nestbc_n, 'n', 's' ) 5228 nzt_topo_nestbc_n, & 5229 kcto, iflo, ifuo, jflo, jfuo, & 5230 kflo, kfuo, ijkfc_s, 'n', 's' ) 4965 5231 ENDDO 4966 5232 ENDIF … … 4971 5237 !-- All PEs are top-border PEs 4972 5238 CALL pmci_interp_tril_t( u, uc, icu, jco, kco, r1xu, r2xu, r1yo, & 4973 r2yo, r1zo, r2zo, 'u' ) 5239 r2yo, r1zo, r2zo, kcto, iflu, ifuu, & 5240 jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' ) 4974 5241 CALL pmci_interp_tril_t( v, vc, ico, jcv, kco, r1xo, r2xo, r1yv, & 4975 r2yv, r1zo, r2zo, 'v' ) 5242 r2yv, r1zo, r2zo, kcto, iflo, ifuo, & 5243 jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' ) 4976 5244 CALL pmci_interp_tril_t( w, wc, ico, jco, kcw, r1xo, r2xo, r1yo, & 4977 r2yo, r1zw, r2zw, 'w' ) 5245 r2yo, r1zw, r2zw, kctw, iflo, ifuo, & 5246 jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' ) 4978 5247 4979 5248 IF ( ( rans_mode_parent .AND. rans_mode ) .OR. & 4980 5249 ( .NOT. rans_mode_parent .AND. .NOT. rans_mode .AND. & 4981 5250 .NOT. constant_diffusion ) ) THEN 4982 CALL pmci_interp_tril_t( e, ec, ico, jco, kco, r1xo, r2xo, r1yo, & 4983 r2yo, r1zo, r2zo, 'e' ) 5251 ! CALL pmci_interp_tril_t( e, ec, ico, jco, kco, r1xo, r2xo, r1yo, & 5252 ! r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5253 ! jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' ) 5254 ! 5255 !-- Interpolation of e is replaced by the Neumann condition. 5256 e(nzt+1,nys:nyn,nxl:nxr) = e(nzt,nys:nyn,nxl:nxr) 5257 4984 5258 ENDIF 4985 5259 4986 5260 IF ( rans_mode_parent .AND. rans_mode .AND. rans_tke_e ) THEN 4987 5261 CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo, & 4988 r1yo, r2yo, r1zo, r2zo, 's' ) 5262 r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5263 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 4989 5264 ENDIF 4990 5265 4991 5266 IF ( .NOT. neutral ) THEN 4992 CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, & 4993 r2yo, r1zo, r2zo, 's' ) 5267 CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, & 5268 r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5269 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 4994 5270 ENDIF 4995 5271 … … 4997 5273 4998 5274 CALL pmci_interp_tril_t( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, & 4999 r2yo, r1zo, r2zo, 's' ) 5275 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5276 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5000 5277 5001 5278 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN 5002 5279 5003 5280 CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,& 5004 r2yo, r1zo, r2zo, 's' ) 5281 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5282 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5005 5283 5006 5284 CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,& 5007 r2yo, r1zo, r2zo, 's' ) 5285 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5286 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5008 5287 5009 5288 ENDIF … … 5013 5292 5014 5293 CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,& 5015 r2yo, r1zo, r2zo, 's' ) 5294 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5295 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5016 5296 5017 5297 CALL pmci_interp_tril_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,& 5018 r2yo, r1zo, r2zo, 's' ) 5298 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5299 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5019 5300 5020 5301 ENDIF … … 5024 5305 IF ( passive_scalar ) THEN 5025 5306 CALL pmci_interp_tril_t( s, sc, ico, jco, kco, r1xo, r2xo, r1yo, & 5026 r2yo, r1zo, r2zo, 's' ) 5307 r2yo, r1zo, r2zo, kcto, iflo, ifuo, & 5308 jflo, jfuo, kflo, kfuo, ijkfc_s, 's' ) 5027 5309 ENDIF 5028 5310 … … 5033 5315 ico, jco, kco, r1xo, r2xo, & 5034 5316 r1yo, r2yo, r1zo, r2zo, & 5035 's' ) 5317 kcto, iflo, ifuo, jflo, jfuo, & 5318 kflo, kfuo, ijkfc_s, 's' ) 5036 5319 ENDDO 5037 5320 ENDIF … … 5050 5333 INTEGER(iwp) :: n !< running index for number of chemical species 5051 5334 5052 5053 5054 CALL pmci_anterp_tophat( u, uc, kctu, iflu, ifuu, jflo, jfuo, kflo, & 5335 CALL pmci_anterp_tophat( u, uc, kcto, iflu, ifuu, jflo, jfuo, kflo, & 5055 5336 kfuo, ijkfc_u, 'u' ) 5056 CALL pmci_anterp_tophat( v, vc, kct u, iflo, ifuo, jflv, jfuv, kflo, &5337 CALL pmci_anterp_tophat( v, vc, kcto, iflo, ifuo, jflv, jfuv, kflo, & 5057 5338 kfuo, ijkfc_v, 'v' ) 5058 5339 CALL pmci_anterp_tophat( w, wc, kctw, iflo, ifuo, jflo, jfuo, kflw, & … … 5062 5343 !-- RANS mode. 5063 5344 IF ( rans_mode_parent .AND. rans_mode ) THEN 5064 CALL pmci_anterp_tophat( e, ec, kct u, iflo, ifuo, jflo, jfuo, kflo, &5345 CALL pmci_anterp_tophat( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, & 5065 5346 kfuo, ijkfc_s, 'e' ) 5066 5347 ! 5067 5348 !-- Anterpolation of dissipation rate only if TKE-e closure is applied. 5068 5349 IF ( rans_tke_e ) THEN 5069 CALL pmci_anterp_tophat( diss, dissc, kct u, iflo, ifuo, jflo, jfuo,&5350 CALL pmci_anterp_tophat( diss, dissc, kcto, iflo, ifuo, jflo, jfuo,& 5070 5351 kflo, kfuo, ijkfc_s, 'diss' ) 5071 5352 ENDIF … … 5074 5355 5075 5356 IF ( .NOT. neutral ) THEN 5076 CALL pmci_anterp_tophat( pt, ptc, kct u, iflo, ifuo, jflo, jfuo, kflo, &5357 CALL pmci_anterp_tophat( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, kflo, & 5077 5358 kfuo, ijkfc_s, 'pt' ) 5078 5359 ENDIF … … 5080 5361 IF ( humidity ) THEN 5081 5362 5082 CALL pmci_anterp_tophat( q, q_c, kct u, iflo, ifuo, jflo, jfuo, kflo, &5363 CALL pmci_anterp_tophat( q, q_c, kcto, iflo, ifuo, jflo, jfuo, kflo, & 5083 5364 kfuo, ijkfc_s, 'q' ) 5084 5365 5085 5366 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN 5086 5367 5087 CALL pmci_anterp_tophat( qc, qcc, kct u, iflo, ifuo, jflo, jfuo, &5368 CALL pmci_anterp_tophat( qc, qcc, kcto, iflo, ifuo, jflo, jfuo, & 5088 5369 kflo, kfuo, ijkfc_s, 'qc' ) 5089 5370 5090 CALL pmci_anterp_tophat( nc, ncc, kct u, iflo, ifuo, jflo, jfuo, &5371 CALL pmci_anterp_tophat( nc, ncc, kcto, iflo, ifuo, jflo, jfuo, & 5091 5372 kflo, kfuo, ijkfc_s, 'nc' ) 5092 5373 … … 5095 5376 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN 5096 5377 5097 CALL pmci_anterp_tophat( qr, qrc, kct u, iflo, ifuo, jflo, jfuo, &5378 CALL pmci_anterp_tophat( qr, qrc, kcto, iflo, ifuo, jflo, jfuo, & 5098 5379 kflo, kfuo, ijkfc_s, 'qr' ) 5099 5380 5100 CALL pmci_anterp_tophat( nr, nrc, kct u, iflo, ifuo, jflo, jfuo, &5381 CALL pmci_anterp_tophat( nr, nrc, kcto, iflo, ifuo, jflo, jfuo, & 5101 5382 kflo, kfuo, ijkfc_s, 'nr' ) 5102 5383 … … 5106 5387 5107 5388 IF ( passive_scalar ) THEN 5108 CALL pmci_anterp_tophat( s, sc, kct u, iflo, ifuo, jflo, jfuo, kflo, &5389 CALL pmci_anterp_tophat( s, sc, kcto, iflo, ifuo, jflo, jfuo, kflo, & 5109 5390 kfuo, ijkfc_s, 's' ) 5110 5391 ENDIF … … 5114 5395 CALL pmci_anterp_tophat( chem_species(n)%conc, & 5115 5396 chem_spec_c(:,:,:,n), & 5116 kct u, iflo, ifuo, jflo, jfuo, kflo, &5397 kcto, iflo, ifuo, jflo, jfuo, kflo, & 5117 5398 kfuo, ijkfc_s, 's' ) 5118 5399 ENDDO … … 5125 5406 SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, & 5126 5407 r2z, logc, logc_ratio, logc_kbounds, & 5127 nzt_topo_nestbc, edge, var ) 5408 nzt_topo_nestbc, & 5409 kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc, & 5410 edge, var ) 5128 5411 ! 5129 5412 !-- Interpolation of ghost-node values used as the child-domain boundary … … 5132 5415 5133 5416 IMPLICIT NONE 5134 5135 INTEGER(iwp) :: nzt_topo_nestbc !<5136 5417 5137 5418 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & … … 5152 5433 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !< 5153 5434 INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nys:nyn), & 5154 INTENT(IN) :: logc !< 5155 INTEGER(iwp), DIMENSION(1:2,nys:nyn), INTENT(IN) :: logc_kbounds !< 5435 INTENT(IN) :: logc !< 5436 INTEGER(iwp), DIMENSION(1:2,nys:nyn), INTENT(IN) :: logc_kbounds !< 5437 5438 INTEGER(iwp) :: kct 5439 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !< Indicates start index of child cells belonging to certain parent cell - x direction 5440 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !< Indicates end index of child cells belonging to certain parent cell - x direction 5441 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !< Indicates start index of child cells belonging to certain parent cell - y direction 5442 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !< Indicates start index of child cells belonging to certain parent cell - y direction 5443 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 5444 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 5445 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 5446 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 5447 5448 !AH INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 5449 INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 5450 5451 INTEGER(iwp) :: nzt_topo_nestbc !< 5156 5452 5157 5453 CHARACTER(LEN=1), INTENT(IN) :: edge !< … … 5159 5455 5160 5456 INTEGER(iwp) :: i !< 5457 INTEGER(iwp) :: ia !< 5161 5458 INTEGER(iwp) :: ib !< 5162 5459 INTEGER(iwp) :: ibgp !< 5460 INTEGER(iwp) :: ijk !< 5461 INTEGER(iwp) :: iw !< 5163 5462 INTEGER(iwp) :: j !< 5164 5463 INTEGER(iwp) :: jco !< 5165 5464 INTEGER(iwp) :: jcorr !< 5166 5465 INTEGER(iwp) :: jinc !< 5466 INTEGER(iwp) :: jw !< 5167 5467 INTEGER(iwp) :: j1 !< 5168 5468 INTEGER(iwp) :: k !< … … 5174 5474 INTEGER(iwp) :: m !< 5175 5475 INTEGER(iwp) :: n !< 5176 5476 INTEGER(iwp) :: kbc !< 5477 INTEGER(iwp) :: var_flag !< 5478 5479 REAL(wp) :: cellsum !< 5480 REAL(wP) :: cellsumd !< 5177 5481 REAL(wp) :: fkj !< 5178 5482 REAL(wp) :: fkjp !< … … 5181 5485 REAL(wp) :: fk !< 5182 5486 REAL(wp) :: fkp !< 5183 5487 REAL(wp) :: rcorr !< 5488 REAL(wp) :: rcorr_ijk !< 5489 5184 5490 ! 5185 5491 !-- Check which edge is to be handled … … 5198 5504 ib = nxr + 2 5199 5505 ENDIF 5200 5506 5507 IF ( var == 'u' ) THEN 5508 var_flag = 1 5509 ELSEIF ( var == 'v' ) THEN 5510 var_flag = 2 5511 ELSEIF ( var == 'w' ) THEN 5512 var_flag = 3 5513 ELSE 5514 var_flag = 0 5515 ENDIF 5516 5201 5517 DO j = nys, nyn+1 5202 5518 ! 5203 5519 !-- Determine vertical index of topography top at grid point (j,i) 5204 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )5205 5206 DO k = k_wall, nzt+15520 !AH k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 5521 5522 DO k = nzb, nzt+1 !k_wall, nzt+1 5207 5523 l = ic(i) 5208 5524 m = jc(j) … … 5235 5551 DO kcorr = 0, ncorr - 1 5236 5552 kco = k + kcorr 5237 f(kco,j,i) = logc_ratio(1,kcorr,k,j) * f(k1,j,i)5553 !AH f(kco,j,i) = logc_ratio(1,kcorr,k,j) * f(k1,j,i) 5238 5554 ENDDO 5239 5555 ENDIF … … 5261 5577 jco = j + jinc * jcorr 5262 5578 IF ( jco >= nys .AND. jco <= nyn ) THEN 5263 f(k,jco,i) = logc_ratio(2,jcorr,k,j) * f(k,j1,i)5579 !AH f(k,jco,i) = logc_ratio(2,jcorr,k,j) * f(k,j1,i) 5264 5580 ENDIF 5265 5581 ENDDO … … 5282 5598 DO kcorr = 0, ncorr-1 5283 5599 kco = k + kcorr 5284 f(kco,jco,i) = 0.5_wp * ( logc_ratio(1,kcorr,k,j) * &5285 f(k1,j,i) &5286 + logc_ratio(2,jcorr,k,j) * &5287 f(k,j1,i) )5600 !AH f(kco,jco,i) = 0.5_wp * ( logc_ratio(1,kcorr,k,j) * & 5601 !AH f(k1,j,i) & 5602 !AH + logc_ratio(2,jcorr,k,j) * & 5603 !AH f(k,j1,i) ) 5288 5604 ENDDO 5289 5605 ENDIF … … 5295 5611 ENDIF ! ( topography /= 'flat' ) 5296 5612 ! 5297 !-- Rescale if f is the TKE. 5298 IF ( var == 'e') THEN 5299 IF ( edge == 'l' ) THEN 5300 DO j = nys, nyn + 1 5301 ! 5302 !-- Determine vertical index of topography top at grid point (j,i) 5303 k_wall = get_topography_top_index_ji( j, i, 's' ) 5304 5305 DO k = k_wall, nzt + 1 5306 f(k,j,i) = tkefactor_l(k,j) * f(k,j,i) 5613 !-- Apply the reversibility correction to the boundary-normal velocity- 5614 !-- component u and the scalars. It must not be applied to the boundary- 5615 !-- tangential velocity components v and w because their 2-D anterpolation 5616 !-- cells do not cover all the child-grid nodes on the boundary. 5617 IF ( .NOT. ( ( var == 'v' ) .OR. ( var == 'w' ) ) ) THEN 5618 l = ic(i) 5619 DO m = jcs, jcn 5620 DO n = 0, kct+1 5621 ijk = 1 5622 cellsum = 0.0_wp 5623 cellsumd = 0.0_wp 5624 ! 5625 !-- Note that the index name i must not be used here as a loop 5626 !-- index name since i is the constant boundary index, hence 5627 !-- the name ia. 5628 DO ia = ifl(l), ifu(l) 5629 DO j = jfl(m), jfu(m) 5630 DO k = kfl(n), kfu(n) 5631 cellsum = cellsum + MERGE( f(k,j,ia), 0.0_wp, & 5632 BTEST( wall_flags_0(k,j,ia), var_flag ) ) 5633 celltmpd(ijk) = ABS( fc(n,m,l) - f(k,j,ia) ) 5634 cellsumd = cellsumd + MERGE( celltmpd(ijk), & 5635 0.0_wp, BTEST( wall_flags_0(k,j,ia), var_flag ) ) 5636 ijk = ijk + 1 5637 ENDDO 5638 ENDDO 5307 5639 ENDDO 5308 ENDDO 5309 ELSEIF ( edge == 'r' ) THEN 5310 DO j = nys, nyn+1 5311 ! 5312 !-- Determine vertical index of topography top at grid point (j,i) 5313 k_wall = get_topography_top_index_ji( j, i, 's' ) 5314 5315 DO k = k_wall, nzt+1 5316 f(k,j,i) = tkefactor_r(k,j) * f(k,j,i) 5640 5641 IF ( ijkfc(n,m,l) /= 0 ) THEN 5642 cellsum = cellsum / REAL( ijkfc(n,m,l), KIND=wp ) 5643 rcorr = fc(n,m,l) - cellsum 5644 cellsumd = cellsumd / REAL( ijkfc(n,m,l), KIND=wp ) 5645 ELSE 5646 cellsum = 0.0_wp 5647 rcorr = 0.0_wp 5648 cellsumd = 1.0_wp 5649 celltmpd = 1.0_wp 5650 ENDIF 5651 ! 5652 !-- Distribute the correction term to the child nodes according to 5653 !-- their relative difference to the parent value such that the 5654 !-- node with the largest difference gets the largest share of the 5655 !-- correction. The distribution is skipped if rcorr is negligibly 5656 !-- small in order to avoid division by zero. 5657 IF ( ABS(rcorr) < 0.000001_wp ) THEN 5658 cellsumd = 1.0_wp 5659 celltmpd = 1.0_wp 5660 ENDIF 5661 5662 ijk = 1 5663 DO ia = ifl(l), ifu(l) 5664 DO j = jfl(m), jfu(m) 5665 DO k = kfl(n), kfu(n) 5666 rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 5667 f(k,j,ia) = f(k,j,ia) + rcorr_ijk 5668 ijk = ijk + 1 5669 ENDDO 5670 ENDDO 5317 5671 ENDDO 5318 ENDDO 5319 ENDIF 5320 ENDIF 5672 5673 ENDDO ! n 5674 ENDDO ! m 5675 5676 ENDIF ! var not v or w 5321 5677 ! 5322 5678 !-- Store the boundary values also into the other redundant ghost node layers. 5323 !-- Please note,in case of only one ghost node layer, e.g. for the PW5679 !-- Note that in case of only one ghost node layer, e.g. for the PW 5324 5680 !-- scheme, the following loops will not be entered. 5325 5681 IF ( edge == 'l' ) THEN … … 5339 5695 SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, & 5340 5696 r2z, logc, logc_ratio, logc_kbounds, & 5341 nzt_topo_nestbc, edge, var ) 5697 nzt_topo_nestbc, & 5698 kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc, & 5699 edge, var ) 5342 5700 5343 5701 ! … … 5347 5705 5348 5706 IMPLICIT NONE 5349 5350 INTEGER(iwp) :: nzt_topo_nestbc !<5351 5707 5352 5708 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & … … 5370 5726 INTEGER(iwp), DIMENSION(1:2,nxl:nxr), INTENT(IN) :: logc_kbounds !< 5371 5727 5728 INTEGER(iwp) :: kct 5729 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !< Indicates start index of child cells belonging to certain parent cell - x direction 5730 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !< Indicates end index of child cells belonging to certain parent cell - x direction 5731 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !< Indicates start index of child cells belonging to certain parent cell - y direction 5732 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !< Indicates start index of child cells belonging to certain parent cell - y direction 5733 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 5734 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 5735 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 5736 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 5737 !AH INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 5738 INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 5739 5740 INTEGER(iwp) :: nzt_topo_nestbc !< 5741 5372 5742 CHARACTER(LEN=1), INTENT(IN) :: edge !< 5373 5743 CHARACTER(LEN=1), INTENT(IN) :: var !< … … 5377 5747 INTEGER(iwp) :: icorr !< 5378 5748 INTEGER(iwp) :: ico !< 5749 INTEGER(iwp) :: ijk !< 5379 5750 INTEGER(iwp) :: i1 !< 5380 5751 INTEGER(iwp) :: j !< 5752 INTEGER(iwp) :: ja !< 5381 5753 INTEGER(iwp) :: jb !< 5382 5754 INTEGER(iwp) :: jbgp !< 5383 5755 INTEGER(iwp) :: k !< 5384 INTEGER(iwp) :: k_wall 5756 INTEGER(iwp) :: k_wall !< vertical index of topography top 5385 5757 INTEGER(iwp) :: kcorr !< 5386 5758 INTEGER(iwp) :: kco !< … … 5389 5761 INTEGER(iwp) :: m !< 5390 5762 INTEGER(iwp) :: n !< 5391 5763 INTEGER(iwp) :: var_flag !< 5764 5765 REAL(wp) :: cellsum !< 5766 REAL(wp) :: cellsumd !< 5392 5767 REAL(wp) :: fk !< 5393 5768 REAL(wp) :: fkj !< … … 5396 5771 REAL(wp) :: fkpjp !< 5397 5772 REAL(wp) :: fkp !< 5773 REAL(wp) :: rcorr !< 5774 REAL(wp) :: rcorr_ijk !< 5398 5775 5399 5776 ! … … 5414 5791 ENDIF 5415 5792 5793 IF ( var == 'u' ) THEN 5794 var_flag = 1 5795 ELSEIF ( var == 'v' ) THEN 5796 var_flag = 2 5797 ELSEIF ( var == 'w' ) THEN 5798 var_flag = 3 5799 ELSE 5800 var_flag = 0 5801 ENDIF 5802 5416 5803 DO i = nxl, nxr+1 5417 5804 ! 5418 5805 !-- Determine vertical index of topography top at grid point (j,i) 5419 k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )5420 5421 DO k = k_wall, nzt+15806 !AH k_wall = get_topography_top_index_ji( j, i, TRIM( var ) ) 5807 5808 DO k = nzb, nzt+1 !AH k_wall, nzt+1 5422 5809 l = ic(i) 5423 5810 m = jc(j) … … 5450 5837 DO kcorr = 0, ncorr-1 5451 5838 kco = k + kcorr 5452 f(kco,j,i) = logc_ratio(1,kcorr,k,i) * f(k1,j,i)5839 !AH f(kco,j,i) = logc_ratio(1,kcorr,k,i) * f(k1,j,i) 5453 5840 ENDDO 5454 5841 ENDIF … … 5477 5864 ico = i + iinc * icorr 5478 5865 IF ( ico >= nxl .AND. ico <= nxr ) THEN 5479 f(k,j,ico) = logc_ratio(2,icorr,k,i) * f(k,j,i1)5866 !AH f(k,j,ico) = logc_ratio(2,icorr,k,i) * f(k,j,i1) 5480 5867 ENDIF 5481 5868 ENDDO … … 5498 5885 DO kcorr = 0, ncorr-1 5499 5886 kco = k + kcorr 5500 f(kco,j,ico) = 0.5_wp * ( logc_ratio(1,kcorr,k,i) * &5501 f(k1,j,i) &5502 + logc_ratio(2,icorr,k,i) * &5503 f(k,j,i1) )5887 !AH f(kco,j,ico) = 0.5_wp * ( logc_ratio(1,kcorr,k,i) * & 5888 !AH f(k1,j,i) & 5889 !AH + logc_ratio(2,icorr,k,i) * & 5890 !AH f(k,j,i1) ) 5504 5891 ENDDO 5505 5892 ENDIF … … 5511 5898 ENDIF ! ( topography /= 'flat' ) 5512 5899 ! 5513 !-- Rescale if f is the TKE. 5514 IF ( var == 'e') THEN 5515 IF ( edge == 's' ) THEN 5516 DO i = nxl, nxr + 1 5517 ! 5518 !-- Determine vertical index of topography top at grid point (j,i) 5519 k_wall = get_topography_top_index_ji( j, i, 's' ) 5520 DO k = k_wall, nzt+1 5521 f(k,j,i) = tkefactor_s(k,i) * f(k,j,i) 5900 !-- Apply the reversibility correction to the boundary-normal velocity- 5901 !-- component v and the scalars. It must not be applied to the boundary- 5902 !-- tangential velocity components u and w because their 2-D anterpolation 5903 !-- cells do not cover all the child-grid nodes on the boundary. 5904 IF ( .NOT. ( ( var == 'u' ) .OR. ( var == 'w' ) ) ) THEN 5905 m = jc(j) 5906 DO l = icl, icr 5907 DO n = 0, kct+1 5908 ijk = 1 5909 cellsum = 0.0_wp 5910 cellsumd = 0.0_wp 5911 DO i = ifl(l), ifu(l) 5912 DO ja = jfl(m), jfu(m) 5913 DO k = kfl(n), kfu(n) 5914 cellsum = cellsum + MERGE( f(k,ja,i), 0.0_wp, & 5915 BTEST( wall_flags_0(k,ja,i), var_flag ) ) 5916 celltmpd(ijk) = ABS( fc(n,m,l) - f(k,ja,i) ) 5917 cellsumd = cellsumd + MERGE( celltmpd(ijk), & 5918 0.0_wp, BTEST( wall_flags_0(k,ja,i), var_flag ) ) 5919 ijk = ijk + 1 5920 ENDDO 5921 ENDDO 5522 5922 ENDDO 5523 ENDDO 5524 ELSEIF ( edge == 'n' ) THEN 5525 DO i = nxl, nxr + 1 5526 ! 5527 !-- Determine vertical index of topography top at grid point (j,i) 5528 k_wall = get_topography_top_index_ji( j, i, 's' ) 5529 DO k = k_wall, nzt+1 5530 f(k,j,i) = tkefactor_n(k,i) * f(k,j,i) 5923 5924 IF ( ijkfc(n,m,l) /= 0 ) THEN 5925 cellsum = cellsum / REAL( ijkfc(n,m,l), KIND=wp ) 5926 rcorr = fc(n,m,l) - cellsum 5927 cellsumd = cellsumd / REAL( ijkfc(n,m,l), KIND=wp ) 5928 ELSE 5929 cellsum = 0.0_wp 5930 rcorr = 0.0_wp 5931 cellsumd = 1.0_wp 5932 celltmpd = 1.0_wp 5933 ENDIF 5934 ! 5935 !-- Distribute the correction term to the child nodes according to 5936 !-- their relative difference to the parent value such that the 5937 !-- node with the largest difference gets the largest share of the 5938 !-- correction. The distribution is skipped if rcorr is negligibly 5939 !-- small in order to avoid division by zero. 5940 IF ( ABS(rcorr) < 0.000001_wp ) THEN 5941 cellsumd = 1.0_wp 5942 celltmpd = 1.0_wp 5943 ENDIF 5944 5945 ijk = 1 5946 DO i = ifl(l), ifu(l) 5947 DO ja = jfl(m), jfu(m) 5948 DO k = kfl(n), kfu(n) 5949 rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 5950 f(k,ja,i) = f(k,ja,i) + rcorr_ijk 5951 ijk = ijk + 1 5952 ENDDO 5953 ENDDO 5531 5954 ENDDO 5532 ENDDO 5533 ENDIF 5534 ENDIF 5955 5956 ENDDO ! n 5957 ENDDO ! l 5958 5959 ENDIF ! var not u or w 5535 5960 ! 5536 5961 !-- Store the boundary values also into the other redundant ghost node layers. 5537 !-- Please note,in case of only one ghost node layer, e.g. for the PW5962 !-- Note that in case of only one ghost node layer, e.g. for the PW 5538 5963 !-- scheme, the following loops will not be entered. 5539 5964 IF ( edge == 's' ) THEN … … 5551 5976 5552 5977 5553 SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, & 5554 r2z, var ) 5978 SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, & 5979 r1z, r2z, kct, ifl, ifu, jfl, jfu, kfl, kfu, & 5980 ijkfc, var ) 5555 5981 5556 5982 ! … … 5576 6002 INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) :: kc !< 5577 6003 6004 INTEGER(iwp) :: kct 6005 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifl !< Indicates start index of child cells belonging to certain parent cell - x direction 6006 INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) :: ifu !< Indicates end index of child cells belonging to certain parent cell - x direction 6007 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !< Indicates start index of child cells belonging to certain parent cell - y direction 6008 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !< Indicates start index of child cells belonging to certain parent cell - y direction 6009 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 6010 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 6011 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 6012 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 6013 !AH INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 6014 INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 6015 5578 6016 CHARACTER(LEN=1), INTENT(IN) :: var !< 5579 6017 … … 5581 6019 INTEGER(iwp) :: ib !< 5582 6020 INTEGER(iwp) :: ie !< 6021 INTEGER(iwp) :: ijk !< 5583 6022 INTEGER(iwp) :: j !< 5584 INTEGER(iwp) :: jb 5585 INTEGER(iwp) :: je 6023 INTEGER(iwp) :: jb !< 6024 INTEGER(iwp) :: je !< 5586 6025 INTEGER(iwp) :: k !< 6026 INTEGER(iwp) :: ka !< 5587 6027 INTEGER(iwp) :: l !< 5588 6028 INTEGER(iwp) :: m !< 5589 6029 INTEGER(iwp) :: n !< 6030 INTEGER(iwp) :: var_flag !< 5590 6031 6032 REAL(wp) :: cellsum !< 6033 REAL(wp) :: cellsumd !< 5591 6034 REAL(wp) :: fk !< 5592 6035 REAL(wp) :: fkj !< … … 5595 6038 REAL(wp) :: fkpjp !< 5596 6039 REAL(wp) :: fkp !< 5597 5598 6040 REAL(wp) :: rcorr !< 6041 REAL(wp) :: rcorr_ijk !< 6042 6043 5599 6044 IF ( var == 'w' ) THEN 5600 6045 k = nzt … … 5610 6055 jb = nys-1 5611 6056 je = nyn+1 5612 ! 5613 !-- The exceedings must not be made past the outer edges in 5614 !-- case of pure vertical nesting. 5615 IF ( nesting_mode == 'vertical' ) THEN 5616 IF ( nxl == 0 ) ib = nxl 5617 IF ( nxr == nx ) ie = nxr 5618 IF ( nys == 0 ) jb = nys 5619 IF ( nyn == ny ) je = nyn 6057 6058 IF ( var == 'u' ) THEN 6059 var_flag = 1 6060 ELSEIF ( var == 'v' ) THEN 6061 var_flag = 2 6062 ELSEIF ( var == 'w' ) THEN 6063 var_flag = 3 6064 ELSE 6065 var_flag = 0 5620 6066 ENDIF 5621 6067 … … 5635 6081 ENDDO 5636 6082 ! 5637 !-- Just fill up the second ghost-node layer for w. 6083 !-- Apply the reversibility correction to the boundary-normal velocity- 6084 !-- component w and scalars. It must not be applied to the boundary- 6085 !-- tangential velocity components u and v because their 2-D anterpolation 6086 !-- cells do not cover all the child-grid nodes on the boundary. 6087 IF ( .NOT. ( ( var == 'u' ) .OR. ( var == 'v' ) ) ) THEN 6088 IF ( var == 'w' ) THEN 6089 n = kc(k) 6090 ELSE 6091 n = kc(k) + 1 6092 ENDIF 6093 6094 DO l = icl, icr 6095 DO m = jcs, jcn 6096 ijk = 1 6097 cellsum = 0.0_wp 6098 cellsumd = 0.0_wp 6099 DO i = ifl(l), ifu(l) 6100 DO j = jfl(m), jfu(m) 6101 DO ka = kfl(n), kfu(n) 6102 cellsum = cellsum + MERGE( f(ka,j,i), 0.0_wp, & 6103 BTEST( wall_flags_0(ka,j,i), var_flag ) ) 6104 celltmpd(ijk) = ABS( fc(n,m,l) - f(ka,j,i) ) 6105 cellsumd = cellsumd + MERGE( celltmpd(ijk), & 6106 0.0_wp, BTEST( wall_flags_0(ka,j,i), var_flag ) ) 6107 ijk = ijk + 1 6108 ENDDO 6109 ENDDO 6110 ENDDO 6111 6112 IF ( ijkfc(n,m,l) /= 0 ) THEN 6113 cellsum = cellsum / REAL( ijkfc(n,m,l), KIND=wp ) 6114 rcorr = fc(n,m,l) - cellsum 6115 cellsumd = cellsumd / REAL( ijkfc(n,m,l), KIND=wp ) 6116 ELSE 6117 cellsum = 0.0_wp 6118 rcorr = 0.0_wp 6119 cellsumd = 1.0_wp 6120 celltmpd = 1.0_wp 6121 ENDIF 6122 6123 IF ( ABS(rcorr) < 0.000001_wp ) THEN 6124 cellsumd = 1.0_wp 6125 celltmpd = 1.0_wp 6126 ENDIF 6127 6128 ijk = 1 6129 DO i = ifl(l), ifu(l) 6130 DO j = jfl(m), jfu(m) 6131 DO ka = kfl(n), kfu(n) 6132 rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 6133 f(ka,j,i) = f(ka,j,i) + rcorr_ijk 6134 ijk = ijk + 1 6135 ENDDO 6136 ENDDO 6137 ENDDO 6138 6139 ENDDO ! m 6140 ENDDO ! l 6141 6142 ENDIF ! var not u or v 6143 ! 6144 !-- Just fill up the redundant second ghost-node layer for w. 5638 6145 IF ( var == 'w' ) THEN 5639 6146 f(nzt+1,:,:) = f(nzt,:,:) 5640 ENDIF5641 !5642 !-- Rescale if f is the TKE.5643 !-- It is assumed that the bottom surface never reaches the top boundary of a5644 !-- nest domain.5645 IF ( var == 'e' ) THEN5646 DO i = nxl, nxr5647 DO j = nys, nyn5648 f(k,j,i) = tkefactor_t(j,i) * f(k,j,i)5649 ENDDO5650 ENDDO5651 6147 ENDIF 5652 6148 … … 5667 6163 CHARACTER(LEN=*), INTENT(IN) :: var !< identifyer for treated variable 5668 6164 6165 !AH INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 6166 INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box 6167 5669 6168 INTEGER(iwp) :: i !< Running index x-direction - fine-grid 6169 INTEGER(iwp) :: icla !< Left boundary index for anterpolation along x 6170 INTEGER(iwp) :: icra !< Right boundary index for anterpolation along x 5670 6171 INTEGER(iwp) :: ii !< Running index x-direction - coarse grid 5671 INTEGER(iwp) :: iclp !< Left boundary index for anterpolation along x5672 INTEGER(iwp) :: icrm !< Right boundary index for anterpolation along x5673 6172 INTEGER(iwp) :: j !< Running index y-direction - fine-grid 6173 INTEGER(iwp) :: jcna !< North boundary index for anterpolation along y 6174 INTEGER(iwp) :: jcsa !< South boundary index for anterpolation along y 5674 6175 INTEGER(iwp) :: jj !< Running index y-direction - coarse grid 5675 INTEGER(iwp) :: jcnm !< North boundary index for anterpolation along y5676 INTEGER(iwp) :: jcsp !< South boundary index for anterpolation along y5677 6176 INTEGER(iwp) :: k !< Running index z-direction - fine-grid 6177 INTEGER(iwp) :: kcb = 0 !< Bottom boundary index for anterpolation along z 5678 6178 INTEGER(iwp) :: kk !< Running index z-direction - coarse grid 5679 INTEGER(iwp) :: kcb = 0 !< Bottom boundary index for anterpolation along z5680 6179 INTEGER(iwp) :: var_flag !< bit number used to flag topography on respective grid 5681 6180 … … 5686 6185 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfl !< Indicates start index of child cells belonging to certain parent cell - y direction 5687 6186 INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) :: jfu !< Indicates start index of child cells belonging to certain parent cell - y direction 5688 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction5689 INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction5690 5691 INTEGER(iwp), DIMENSION(0: kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box6187 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 6188 !AH INTEGER(iwp), DIMENSION(0:kct), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 6189 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfl !< Indicates start index of child cells belonging to certain parent cell - z direction 6190 INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN) :: kfu !< Indicates start index of child cells belonging to certain parent cell - z direction 5692 6191 5693 6192 REAL(wp) :: cellsum !< sum of respective child cells belonging to parent cell … … 5698 6197 5699 6198 ! 5700 !-- Initialize the index bounds for anterpolation 5701 iclp = icl 5702 icrm = icr 5703 jcsp = jcs 5704 jcnm = jcn 5705 kcb = 0 5706 ! 5707 !-- Define the index bounds iclp, icrm, jcsp and jcnm. 6199 !-- Define the index bounds icla, icra, jcsa and jcna. 5708 6200 !-- Note that kcb is simply zero and kct enters here as a parameter and it is 5709 6201 !-- determined in pmci_init_anterp_tophat. … … 5714 6206 !-- child domain, leading to increased velocity variances. A more 5715 6207 !-- comprehensive explanation for this is still pending. 5716 IF ( nesting_mode == 'vertical' ) THEN 5717 IF ( bc_dirichlet_l ) THEN 5718 iclp = icl + nhll 5719 ENDIF 5720 IF ( bc_dirichlet_r ) THEN 5721 icrm = icr - nhlr 5722 ENDIF 5723 IF ( bc_dirichlet_s ) THEN 5724 jcsp = jcs + nhls 5725 ENDIF 5726 IF ( bc_dirichlet_n ) THEN 5727 jcnm = jcn - nhln 5728 ENDIF 5729 ELSE 6208 icla = coarse_bound_anterp(1) 6209 icra = coarse_bound_anterp(2) 6210 jcsa = coarse_bound_anterp(3) 6211 jcna = coarse_bound_anterp(4) 6212 kcb = 0 6213 IF ( nesting_mode /= 'vertical' ) THEN 5730 6214 IF ( bc_dirichlet_l ) THEN 5731 6215 IF ( var == 'u' ) THEN 5732 icl p = icl + nhll + 1 + 16216 icla = coarse_bound_anterp(1) + 2 5733 6217 ELSE 5734 icl p = icl + nhll+ 16218 icla = coarse_bound_anterp(1) + 1 5735 6219 ENDIF 5736 6220 ENDIF 5737 6221 IF ( bc_dirichlet_r ) THEN 5738 icr m = icr - nhlr- 16222 icra = coarse_bound_anterp(2) - 1 5739 6223 ENDIF 5740 6224 5741 6225 IF ( bc_dirichlet_s ) THEN 5742 6226 IF ( var == 'v' ) THEN 5743 jcs p = jcs + nhls + 1 + 16227 jcsa = coarse_bound_anterp(3) + 2 5744 6228 ELSE 5745 jcs p = jcs + nhls+ 16229 jcsa = coarse_bound_anterp(3) + 1 5746 6230 ENDIF 5747 6231 ENDIF 5748 6232 IF ( bc_dirichlet_n ) THEN 5749 jcn m = jcn - nhln- 16233 jcna = coarse_bound_anterp(4) - 1 5750 6234 ENDIF 5751 6235 ENDIF 6236 6237 ! write(9,"('pmci_anterp_tophat: ',4(e12.5,2x))") & 6238 ! cg%coord_x(icla), cg%coord_y(jcsa), cg%coord_x(icra), cg%coord_y(jcna) 6239 ! flush(9) 5752 6240 ! 5753 6241 !-- Set masking bit for topography flags … … 5761 6249 var_flag = 0 5762 6250 ENDIF 5763 5764 6251 ! 5765 6252 !-- Note that ii, jj, and kk are coarse-grid indices and i,j, and k 5766 6253 !-- are fine-grid indices. 5767 DO ii = icl p, icrm5768 DO jj = jcs p, jcnm6254 DO ii = icla, icra 6255 DO jj = jcsa, jcna 5769 6256 ! 5770 6257 !-- For simplicity anterpolate within buildings and under elevated 5771 6258 !-- terrain too 5772 DO kk = kcb, kct - 1 5773 6259 DO kk = kcb, kct - 1 5774 6260 cellsum = 0.0_wp 5775 6261 DO i = ifl(ii), ifu(ii) … … 5783 6269 ! 5784 6270 !-- Spatial under-relaxation. 5785 fra = frax(ii) * fray(jj) * fraz(kk) 6271 !-- The relaxation buffer zones are no longer needed with 6272 !-- the new reversible interpolation algorithm. 23.19.2018. 6273 ! fra = frax(ii) * fray(jj) * fraz(kk) 5786 6274 ! 5787 6275 !-- In case all child grid points are inside topography, i.e. … … 5791 6279 !-- zero, keep the parent solution at this point. 5792 6280 IF ( ijkfc(kk,jj,ii) /= 0 ) THEN 5793 fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) + & 5794 fra * cellsum / & 5795 REAL( ijkfc(kk,jj,ii), KIND=wp ) 5796 ENDIF 6281 ! fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) + & 6282 ! fra * cellsum / & 6283 ! REAL( ijkfc(kk,jj,ii), KIND=wp ) 6284 fc(kk,jj,ii) = cellsum / REAL( ijkfc(kk,jj,ii), KIND=wp ) 6285 ENDIF 5797 6286 5798 6287 ENDDO -
palm/trunk/SOURCE/time_integration.f90
r3473 r3484 25 25 ! ----------------- 26 26 ! $Id$ 27 ! pmci_ensure_nest_mass_conservation is premanently removed 28 ! 29 ! 3473 2018-10-30 20:50:15Z suehring 27 30 ! new module for virtual measurements introduced 28 31 ! … … 517 520 USE pmc_interface, & 518 521 ONLY: nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans, & 519 pmci_ ensure_nest_mass_conservation, pmci_synchronize522 pmci_synchronize 520 523 521 524 USE progress_bar, & … … 978 981 !-- Set boundary conditions again after interpolation and anterpolation. 979 982 CALL pmci_boundary_conds 980 !981 !-- Correct the w top-BC in nest domains to ensure mass conservation.982 !-- This action must never be done for the root domain. Vertical983 !-- Commented out April 18, 2018 as seemingly unnecessary.984 !-- Will later be completely removed.985 !-- IF ( child_domain ) THEN986 !-- CALL pmci_ensure_nest_mass_conservation987 !-- ENDIF988 989 983 990 984 CALL cpu_log( log_point(60), 'nesting', 'stop' )
Note: See TracChangeset
for help on using the changeset viewer.
