Changeset 2997 for palm/trunk


Ignore:
Timestamp:
Apr 19, 2018 1:35:17 PM (7 years ago)
Author:
suehring
Message:

Mask topography grid points in anterpolation.

File:
1 edited

Legend:

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

    r2984 r2997  
    2525! -----------------
    2626! $Id$
     27! Mask topography grid points in anterpolation
     28!
     29! 2984 2018-04-18 11:51:30Z hellstea
    2730! Bugfix in the log-law correction initialization. Pivot node cannot any more be
    2831! selected from outside the subdomain in order to prevent array under/overflows.
     
    511514!
    512515!-- Child-array indices to be precomputed and stored for anterpolation.
    513     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !<
    514     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !<
    515     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !<
    516     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !<
    517     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !<
    518     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !<
    519     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !<
    520     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !<
    521     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !<
    522     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !<
    523     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !<
    524     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !<
     516    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !< child index indicating left bound of parent grid box on u-grid
     517    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !< child index indicating right bound of parent grid box on u-grid
     518    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !< child index indicating left bound of parent grid box on scalar-grid
     519    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !< child index indicating right bound of parent grid box on scalar-grid
     520    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !< child index indicating south bound of parent grid box on v-grid
     521    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !< child index indicating north bound of parent grid box on v-grid
     522    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !< child index indicating south bound of parent grid box on scalar-grid
     523    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !< child index indicating north bound of parent grid box on scalar-grid
     524    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !< child index indicating lower bound of parent grid box on w-grid
     525    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !< child index indicating upper bound of parent grid box on w-grid
     526    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !< child index indicating lower bound of parent grid box on scalar-grid
     527    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !< child index indicating upper bound of parent grid box on scalar-grid
    525528!
    526529!-- Number of fine-grid nodes inside coarse-grid ij-faces
    527530!-- to be precomputed for anterpolation.
    528     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_u        !<
    529     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_v        !<
    530     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_s        !<
    531     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_w         !<
    532     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_s         !<
     531    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid boxes contribution to a parent grid box, u-grid
     532    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid boxes contribution to a parent grid box, v-grid
     533    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid boxes contribution to a parent grid box, w-grid
     534    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid boxes contribution to a parent grid box, scalar-grid
    533535   
    534536    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !<
     
    29782980
    29792981       INTEGER(iwp) ::  i        !< Fine-grid index
    2980        INTEGER(iwp) ::  ifc_o    !<
    2981        INTEGER(iwp) ::  ifc_u    !<
    29822982       INTEGER(iwp) ::  ii       !< Coarse-grid index
    29832983       INTEGER(iwp) ::  istart   !<
     
    30393039       ALLOCATE( kfuo(0:kctu) )
    30403040
    3041        ALLOCATE( ijfc_u(jcs:jcn,icl:icr) )
    3042        ALLOCATE( ijfc_v(jcs:jcn,icl:icr) )
    3043        ALLOCATE( ijfc_s(jcs:jcn,icl:icr) )
    3044        ALLOCATE( kfc_w(0:kctw) )
    3045        ALLOCATE( kfc_s(0:kctu) )
     3041       ALLOCATE( ijkfc_u(0:kctu,jcs:jcn,icl:icr) )
     3042       ALLOCATE( ijkfc_v(0:kctu,jcs:jcn,icl:icr) )
     3043       ALLOCATE( ijkfc_w(0:kctw,jcs:jcn,icl:icr) )
     3044       ALLOCATE( ijkfc_s(0:kctu,jcs:jcn,icl:icr) )
     3045       ijkfc_u = 0
     3046       ijkfc_v = 0
     3047       ijkfc_w = 0
     3048       ijkfc_s = 0
    30463049!
    30473050!--    i-indices of u for each ii-index value
     
    31713174!--    This information is needed in anterpolation.
    31723175       DO  ii = icl, icr
    3173           ifc_u = ifuu(ii) - iflu(ii) + 1
    3174           ifc_o = ifuo(ii) - iflo(ii) + 1
    31753176          DO  jj = jcs, jcn
    3176              ijfc_u(jj,ii) = ifc_u * ( jfuo(jj) - jflo(jj) + 1 )
    3177              ijfc_v(jj,ii) = ifc_o * ( jfuv(jj) - jflv(jj) + 1 )
    3178              ijfc_s(jj,ii) = ifc_o * ( jfuo(jj) - jflo(jj) + 1 )             
     3177             DO kk = 0, kctu
     3178!
     3179!--             u-component
     3180                DO  i = iflu(ii), ifuu(ii)
     3181                   DO  j = jflo(jj), jfuo(jj)
     3182                      DO  k = kflo(kk), kfuo(kk)
     3183                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii) + MERGE( 1, 0,  &
     3184                                            BTEST( wall_flags_0(k,j,i), 1 ) )
     3185                      ENDDO
     3186                   ENDDO
     3187                ENDDO
     3188!
     3189!--             v-component
     3190                DO  i = iflo(ii), ifuo(ii)
     3191                   DO  j = jflv(jj), jfuv(jj)
     3192                      DO  k = kflo(kk), kfuo(kk)
     3193                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii) + MERGE( 1, 0,  &
     3194                                            BTEST( wall_flags_0(k,j,i), 2 ) )
     3195                      ENDDO
     3196                   ENDDO
     3197                ENDDO
     3198!
     3199!--             scalars
     3200                DO  i = iflo(ii), ifuo(ii)
     3201                   DO  j = jflo(jj), jfuo(jj)
     3202                      DO  k = kflo(kk), kfuo(kk)
     3203                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii) + MERGE( 1, 0,  &
     3204                                            BTEST( wall_flags_0(k,j,i), 0 ) )
     3205                      ENDDO
     3206                   ENDDO
     3207                ENDDO
     3208             ENDDO
     3209
     3210             DO kk = 0, kctw
     3211!
     3212!--             w-component
     3213                DO  i = iflo(ii), ifuo(ii)
     3214                   DO  j = jflo(jj), jfuo(jj)
     3215                      DO  k = kflw(kk), kfuw(kk)
     3216                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii) + MERGE( 1, 0,  &
     3217                                            BTEST( wall_flags_0(k,j,i), 3 ) )
     3218                      ENDDO
     3219                   ENDDO
     3220                ENDDO
     3221             ENDDO
     3222       
    31793223          ENDDO
    3180        ENDDO
    3181        DO kk = 0, kctw
    3182           kfc_w(kk) = kfuw(kk) - kflw(kk) + 1
    3183        ENDDO
    3184        DO kk = 0, kctu
    3185           kfc_s(kk) = kfuo(kk) - kflo(kk) + 1
    31863224       ENDDO
    31873225!
     
    50125050
    50135051      CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
    5014                                kfuo, ijfc_u, kfc_s, 'u' )
     5052                               kfuo, ijkfc_u, 'u' )
    50155053      CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
    5016                                kfuo, ijfc_v, kfc_s, 'v' )
     5054                               kfuo, ijkfc_v, 'v' )
    50175055      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
    5018                                kfuw, ijfc_s, kfc_w, 'w' )
     5056                               kfuw, ijkfc_w, 'w' )
    50195057!
    50205058!--   Anterpolation of TKE and dissipation rate if parent and child are in
     
    50225060      IF ( rans_mode_parent  .AND.  rans_mode )  THEN
    50235061         CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
    5024                                   kfuo, ijfc_s, kfc_s, 'e' )
     5062                                  kfuo, ijkfc_s, 'e' )
    50255063!
    50265064!--      Anterpolation of dissipation rate only if TKE-e closure is applied.
    50275065         IF ( rans_tke_e )  THEN
    50285066            CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
    5029                                      kflo, kfuo, ijfc_s, kfc_s, 'diss' )
     5067                                     kflo, kfuo, ijkfc_s, 'diss' )
    50305068         ENDIF
    50315069
     
    50345072      IF ( .NOT. neutral )  THEN
    50355073         CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo, &
    5036                                   kfuo, ijfc_s, kfc_s, 'pt' )
     5074                                  kfuo, ijkfc_s, 'pt' )
    50375075      ENDIF
    50385076
     
    50405078
    50415079         CALL pmci_anterp_tophat( q, q_c, kctu, iflo, ifuo, jflo, jfuo, kflo,  &
    5042                                   kfuo, ijfc_s, kfc_s, 'q' )
     5080                                  kfuo, ijkfc_s, 'q' )
    50435081
    50445082         IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
    50455083
    50465084            CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
    5047                                      kflo, kfuo, ijfc_s, kfc_s, 'qc' )
     5085                                     kflo, kfuo, ijkfc_s, 'qc' )
    50485086
    50495087            CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
    5050                                      kflo, kfuo, ijfc_s, kfc_s, 'nc' )
     5088                                     kflo, kfuo, ijkfc_s, 'nc' )
    50515089
    50525090         ENDIF
     
    50555093
    50565094            CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
    5057                                      kflo, kfuo, ijfc_s, kfc_s, 'qr' )
     5095                                     kflo, kfuo, ijkfc_s, 'qr' )
    50585096
    50595097            CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &
    5060                                      kflo, kfuo, ijfc_s, kfc_s, 'nr' )
     5098                                     kflo, kfuo, ijkfc_s, 'nr' )
    50615099
    50625100         ENDIF
     
    50665104      IF ( passive_scalar )  THEN
    50675105         CALL pmci_anterp_tophat( s, sc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
    5068                                   kfuo, ijfc_s, kfc_s, 's' )
     5106                                  kfuo, ijkfc_s, 's' )
    50695107      ENDIF
    50705108
     
    50745112                                     chem_spec_c(:,:,:,n),                     &
    50755113                                     kctu, iflo, ifuo, jflo, jfuo, kflo,       &
    5076                                      kfuo, ijfc_s, kfc_s, 's' )
     5114                                     kfuo, ijkfc_s, 's' )
    50775115         ENDDO
    50785116      ENDIF
     
    56255663
    56265664    SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu,   &
    5627                                    ijfc, kfc, var )
     5665                                   ijkfc, var )
    56285666!
    56295667!--    Anterpolation of internal-node values to be used as the parent-domain
     
    56355673
    56365674       CHARACTER(LEN=*), INTENT(IN) ::  var   !< identifyer for treated variable
     5675
     5676       INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
    56375677
    56385678       INTEGER(iwp) ::  i         !< Running index x-direction - fine-grid
     
    56475687       INTEGER(iwp) ::  kk        !< Running index z-direction - coarse grid
    56485688       INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
    5649        INTEGER(iwp) ::  m         !< Running index surface type
    5650        INTEGER(iwp) ::  nfc       !<
     5689       INTEGER(iwp) ::  var_flag  !< bit number used to flag topography on respective grid
    56515690
    56525691       INTEGER(iwp), INTENT(IN) ::  kct   !< Top boundary index for anterpolation along z
    56535692
    5654        INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl         !< Indicates start index of child cells belonging to certain parent cell - x direction
    5655        INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu         !< Indicates end index of child cells belonging to certain parent cell - x direction
    5656        INTEGER(iwp), DIMENSION(jcs:jcn,icl:icr), INTENT(IN) :: ijfc !<
    5657        INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl         !< Indicates start index of child cells belonging to certain parent cell - y direction
    5658        INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu         !< Indicates start index of child cells belonging to certain parent cell - y direction
    5659        INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfc         !<
    5660        INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl         !< Indicates start index of child cells belonging to certain parent cell - z direction
    5661        INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu         !< Indicates start index of child cells belonging to certain parent cell - z direction
     5693       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
     5694       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
     5695       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
     5696       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
     5697       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
     5698       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
    56625699
    56635700       REAL(wp) ::  cellsum   !< sum of respective child cells belonging to parent cell
     
    57195756             jcnm = jcn - nhln - 1
    57205757          ENDIF
    5721        ENDIF     
     5758       ENDIF
     5759!
     5760!--    Set masking bit for topography flags
     5761       IF ( var == 'u' )  THEN
     5762          var_flag = 1
     5763       ELSEIF ( var == 'v' )  THEN
     5764          var_flag = 2
     5765       ELSEIF ( var == 'w' )  THEN
     5766          var_flag = 3
     5767       ELSE
     5768          var_flag = 0
     5769       ENDIF 
     5770
    57225771!
    57235772!--    Note that ii, jj, and kk are coarse-grid indices and i,j, and k
     
    57295778!--          terrain too
    57305779             DO  kk = kcb, kct - 1
    5731 !
    5732 !--             ijfc and kfc are precomputed in pmci_init_anterp_tophat
    5733                 nfc =  ijfc(jj,ii) * kfc(kk)
     5780
    57345781                cellsum = 0.0_wp
    57355782                DO  i = ifl(ii), ifu(ii)
    57365783                   DO  j = jfl(jj), jfu(jj)
    57375784                      DO  k = kfl(kk), kfu(kk)
    5738                          cellsum = cellsum + f(k,j,i)
     5785                         cellsum = cellsum + MERGE( f(k,j,i), 0.0_wp,          &
     5786                                        BTEST( wall_flags_0(k,j,i), var_flag ) )
    57395787                      ENDDO
    57405788                   ENDDO
     
    57435791!--             Spatial under-relaxation.
    57445792                fra  = frax(ii) * fray(jj) * fraz(kk)
    5745                 fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +               &
    5746                                fra * cellsum / REAL( nfc, KIND = wp )
     5793!
     5794!--             In case all child grid points are inside topography, i.e.
     5795!--             cellsum is zero, also parent solution would have zero values at
     5796!--             that point, which may cause problems in particular for the
     5797!--             temperature. Therefore, in case cellsum is zero, keep the
     5798!--             parent solution at this point. 
     5799                fc(kk,jj,ii) = MERGE( ( 1.0_wp - fra ) * fc(kk,jj,ii) +        &
     5800                                      fra * cellsum  /                         &
     5801                                      REAL( ijkfc(kk,jj,ii), KIND = wp ),      &
     5802                                      fc(kk,jj,ii),                            &
     5803                                      cellsum /= 0.0_wp )       
    57475804
    57485805             ENDDO
Note: See TracChangeset for help on using the changeset viewer.