Changeset 2019 for palm


Ignore:
Timestamp:
Sep 30, 2016 1:40:09 PM (8 years ago)
Author:
hellstea
Message:

bugfixes in pmc_interface_mod.f90

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/pmc_interface_mod.f90

    r2005 r2019  
    2121! Current revisions:
    2222! ------------------
    23 !
     23! Bugfixes mainly in determining the anterpolation index bounds. These errors
     24! were detected when first time tested using 3:1 grid-spacing.
    2425!
    2526! Former revisions:
     
    126127   USE arrays_3d,                                                               &
    127128        ONLY:  dzu, dzw, e, e_p, e_1, e_2, pt, pt_p, pt_1, pt_2, q, q_p, q_1,   &
    128                q_2, s, s_2, u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2,  &
    129                zu, zw, z0
     129               q_2, s, s_2, u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1,    &
     130               w_2, zu, zw, z0
    130131#endif
    131132
     
    10901091!--    This is needed in the anterpolation phase
    10911092       nhlr = CEILING( xexr / cg%dx )
    1092        xce  = coord_x(nxr) + xexr
     1093       xce  = coord_x(nxr+1) + xexr
     1094!--    One "extra" layer is taken behind the right boundary
     1095!--    because it may be needed in cases of non-integer grid-spacing ratio
    10931096       DO  i = cg%nx, 0 , -1
    10941097          IF ( cg%coord_x(i) < xce )  THEN
     
    11181121!--    This is needed in the anterpolation phase
    11191122       nhln = CEILING( yexn / cg%dy )
    1120        yce  = coord_y(nyn) + yexn
     1123       yce  = coord_y(nyn+1) + yexn
     1124!--    One "extra" layer is taken behind the north boundary
     1125!--    because it may be needed in cases of non-integer grid-spacing ratio
    11211126       DO  j = cg%ny, 0, -1
    11221127          IF ( cg%coord_y(j) < yce )  THEN
     
    13581363
    13591364!
    1360 !--    First horizontal walls
     1365!--    First horizontal walls. Note that also logc_w_? and logc_ratio_w_? need to
     1366!--    be allocated and initialized here.
    13611367!--    Left boundary
    13621368       IF ( nest_bound_l )  THEN
     
    13641370          ALLOCATE( logc_u_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
    13651371          ALLOCATE( logc_v_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
     1372          ALLOCATE( logc_w_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
    13661373          ALLOCATE( logc_ratio_u_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
    13671374          ALLOCATE( logc_ratio_v_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
     1375          ALLOCATE( logc_ratio_w_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
    13681376          logc_u_l       = 0
    13691377          logc_v_l       = 0
     1378          logc_w_l       = 0
    13701379          logc_ratio_u_l = 1.0_wp
    13711380          logc_ratio_v_l = 1.0_wp
     1381          logc_ratio_w_l = 1.0_wp
    13721382          direction      = 1
    13731383          inc            = 1
     
    14041414!--    Right boundary
    14051415       IF ( nest_bound_r )  THEN
    1406 
     1416           
    14071417          ALLOCATE( logc_u_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
    14081418          ALLOCATE( logc_v_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
     1419          ALLOCATE( logc_w_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
    14091420          ALLOCATE( logc_ratio_u_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
    14101421          ALLOCATE( logc_ratio_v_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
     1422          ALLOCATE( logc_ratio_w_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
    14111423          logc_u_r       = 0
    14121424          logc_v_r       = 0
     1425          logc_w_r       = 0
    14131426          logc_ratio_u_r = 1.0_wp
    14141427          logc_ratio_v_r = 1.0_wp
     1428          logc_ratio_w_r = 1.0_wp
    14151429          direction      = 1
    14161430          inc            = 1
     1431
    14171432          DO  j = nys, nyn
    14181433!
     
    14491464          ALLOCATE( logc_u_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
    14501465          ALLOCATE( logc_v_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
     1466          ALLOCATE( logc_w_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
    14511467          ALLOCATE( logc_ratio_u_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
    14521468          ALLOCATE( logc_ratio_v_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
     1469          ALLOCATE( logc_ratio_w_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
    14531470          logc_u_s       = 0
    14541471          logc_v_s       = 0
     1472          logc_w_s       = 0
    14551473          logc_ratio_u_s = 1.0_wp
    14561474          logc_ratio_v_s = 1.0_wp
     1475          logc_ratio_w_s = 1.0_wp
    14571476          direction      = 1
    14581477          inc            = 1
     
    14921511          ALLOCATE( logc_u_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
    14931512          ALLOCATE( logc_v_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
     1513          ALLOCATE( logc_w_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
    14941514          ALLOCATE( logc_ratio_u_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
    14951515          ALLOCATE( logc_ratio_v_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
     1516          ALLOCATE( logc_ratio_w_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
    14961517          logc_u_n       = 0
    14971518          logc_v_n       = 0
     1519          logc_w_n       = 0
    14981520          logc_ratio_u_n = 1.0_wp
    14991521          logc_ratio_v_n = 1.0_wp
     1522          logc_ratio_w_n = 1.0_wp
    15001523          direction      = 1
    15011524          inc            = 1
     
    15381561          IF ( nest_bound_l )  THEN
    15391562
    1540              ALLOCATE( logc_w_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
    1541              ALLOCATE( logc_ratio_w_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,      &
    1542                                       nys:nyn) )
    1543 
    1544              logc_w_l       = 0
    1545              logc_ratio_w_l = 1.0_wp
    1546              direction      = 2
     1563             direction  = 2
     1564
    15471565             DO  j = nys, nyn
    15481566                DO  k = nzb, nzt_topo_nestbc_l
     
    16251643          IF ( nest_bound_r )  THEN
    16261644
    1627              ALLOCATE( logc_w_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
    1628              ALLOCATE( logc_ratio_w_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,      &
    1629                                       nys:nyn) )
    1630              logc_w_r       = 0
    1631              logc_ratio_w_r = 1.0_wp
    1632              direction      = 2
     1645             direction  = 2
    16331646             i  = nxr + 1
    16341647
     
    17091722          IF ( nest_bound_s )  THEN
    17101723
    1711              ALLOCATE( logc_w_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
    1712              ALLOCATE( logc_ratio_w_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,      &
    1713                                       nxl:nxr) )
    1714              logc_w_s       = 0
    1715              logc_ratio_w_s = 1.0_wp
    1716              direction      = 3
     1724             direction  = 3
    17171725
    17181726             DO  i = nxl, nxr
     
    17951803          IF ( nest_bound_n )  THEN
    17961804
    1797              ALLOCATE( logc_w_n(1:2,nzb:nzt_topo_nestbc_n, nxl:nxr) )
    1798              ALLOCATE( logc_ratio_w_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,      &
    1799                                       nxl:nxr) )
    1800              logc_w_n       = 0
    1801              logc_ratio_w_n = 1.0_wp
    1802              direction      = 3
     1805             direction  = 3
    18031806             j  = nyn + 1
    18041807
     
    21522155!--    index k
    21532156       kk = 0
    2154        DO  WHILE ( cg%zu(kk) < zu(nzt) )
     2157       DO  WHILE ( cg%zu(kk) <= zu(nzt) )
    21552158          kk = kk + 1
    21562159       ENDDO
     
    21582161
    21592162       kk = 0
    2160        DO  WHILE ( cg%zw(kk) < zw(nzt-1) )
     2163       DO  WHILE ( cg%zw(kk) <= zw(nzt-1) )
    21612164          kk = kk + 1
    21622165       ENDDO
     
    21822185!
    21832186!--    i-indices of u for each ii-index value
     2187!--    ii=icr is redundant for anterpolation
    21842188       istart = nxlg
    2185        DO  ii = icl, icr
     2189       DO  ii = icl, icr-1
    21862190          i = istart
    21872191          DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.   &
     
    21902194          ENDDO
    21912195          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
    2192           DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.   &
    2193                       ( i < nxrg ) )
     2196          DO  WHILE ( ( coord_x(i) <= cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.  &
     2197                      ( i < nxrg+1 ) )
    21942198             i = i + 1
    21952199          ENDDO
    2196           ifuu(ii) = MIN( MAX( i, nxlg ), nxrg )
     2200          ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg )
    21972201          istart = iflu(ii)
    21982202       ENDDO
     2203       iflu(icr) = nxrg
     2204       ifuu(icr) = nxrg
    21992205
    22002206!
    22012207!--    i-indices of others for each ii-index value
     2208!--    ii=icr is redundant for anterpolation
    22022209       istart = nxlg
    2203        DO  ii = icl, icr
     2210       DO  ii = icl, icr-1
    22042211          i = istart
    22052212          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) )  .AND.      &
     
    22082215          ENDDO
    22092216          iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
    2210           DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) + cg%dx )     &
    2211                       .AND.  ( i < nxrg ) )
     2217          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx <= cg%coord_x(ii) + cg%dx )    &
     2218                      .AND.  ( i < nxrg+1 ) )
    22122219             i = i + 1
    22132220          ENDDO
    2214           ifuo(ii) = MIN(MAX(i,nxlg),nxrg)
     2221          ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg )
    22152222          istart = iflo(ii)
    22162223       ENDDO
     2224       iflo(icr) = nxrg
     2225       ifuo(icr) = nxrg
    22172226
    22182227!
    22192228!--    j-indices of v for each jj-index value
     2229!--    jj=jcn is redundant for anterpolation
    22202230       jstart = nysg
    2221        DO  jj = jcs, jcn
     2231       DO  jj = jcs, jcn-1
    22222232          j = jstart
    22232233          DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.   &
     
    22262236          ENDDO
    22272237          jflv(jj) = MIN( MAX( j, nysg ), nyng )
    2228           DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.   &
    2229                       ( j < nyng ) )
     2238          DO  WHILE ( ( coord_y(j) <= cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.  &
     2239                      ( j < nyng+1 ) )
    22302240             j = j + 1
    22312241          ENDDO
    2232           jfuv(jj) = MIN( MAX( j, nysg ), nyng )
     2242          jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng )
    22332243          jstart = jflv(jj)
    22342244       ENDDO
    2235 
     2245       jflv(jcn) = nyng
     2246       jfuv(jcn) = nyng
    22362247!
    22372248!--    j-indices of others for each jj-index value
     2249!--    jj=jcn is redundant for anterpolation
    22382250       jstart = nysg
    2239        DO  jj = jcs, jcn
     2251       DO  jj = jcs, jcn-1
    22402252          j = jstart
    22412253          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) )  .AND.      &
     
    22442256          ENDDO
    22452257          jflo(jj) = MIN( MAX( j, nysg ), nyng )
    2246           DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) + cg%dy )     &
    2247                       .AND.  ( j < nyng ) )
     2258          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy <= cg%coord_y(jj) + cg%dy )    &
     2259                      .AND.  ( j < nyng+1 ) )
    22482260             j = j + 1
    22492261          ENDDO
    2250           jfuo(jj) = MIN( MAX( j, nysg ), nyng )
    2251           jstart = jflv(jj)
     2262          jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng )
     2263          jstart = jflo(jj)
    22522264       ENDDO
     2265       jflo(jcn) = nyng
     2266       jfuo(jcn) = nyng
    22532267
    22542268!
     
    22592273       DO  kk = 1, kctw
    22602274          k = kstart
    2261           DO  WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) )  .AND.        &
    2262                       ( k < nzt ) )
     2275          DO  WHILE ( ( zw(k) < cg%zu(kk) )  .AND.  ( k < nzt ) )
    22632276             k = k + 1
    22642277          ENDDO
    22652278          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
    2266           DO  WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) )  .AND.      &
    2267                       ( k < nzt ) )
     2279          DO  WHILE ( ( zw(k) <= cg%zu(kk+1) )  .AND.  ( k < nzt+1 ) )
    22682280             k = k + 1
    22692281          ENDDO
    2270           kfuw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
     2282          kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 )
    22712283          kstart = kflw(kk)
    22722284       ENDDO
     
    22792291       DO  kk = 1, kctu
    22802292          k = kstart
    2281           DO  WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) )  .AND.        &
    2282                       ( k < nzt ) )
     2293          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k < nzt ) )
    22832294             k = k + 1
    22842295          ENDDO
    22852296          kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
    2286           DO  WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) )  .AND.      &
    2287                       ( k < nzt ) )
     2297          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
    22882298             k = k + 1
    22892299          ENDDO
    2290           kfuo(kk) = MIN( MAX( k-1, 1 ), nzt + 1 )
     2300          kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 )
    22912301          kstart = kflo(kk)
    22922302       ENDDO
Note: See TracChangeset for help on using the changeset viewer.