Ignore:
Timestamp:
Aug 16, 2019 1:50:17 PM (2 years ago)
Author:
suehring
Message:

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

File:
1 edited

Legend:

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

    r4167 r4168  
    2828! -----------------
    2929! $Id$
    30 ! Changed behaviour of masked output over surface to follow terrain and ignore
    31 ! buildings (J.Resler, T.Gronemeier)
     30! Replace function get_topography_top_index by topo_top_ind
    3231!
    3332! 4157 2019-08-14 09:19:12Z suehring
     
    691690    USE indices,                                                               &
    692691        ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
    693                nzb, nzt
     692               nzb, nzt, topo_top_ind
    694693
    695694    USE, INTRINSIC :: iso_c_binding
     
    737736
    738737    USE surface_mod,                                                           &
    739         ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
    740                ind_pav_green, ind_veg_wall, ind_wat_win,                       &
     738        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
    741739               surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
    742740               vertical_surfaces_exist
     
    745743
    746744    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
    747 
    748     REAL(wp), PARAMETER ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
    749745
    750746!
     
    35833579!
    35843580!--       Determine minimum topography top index.
    3585           k_topo_l = MINVAL( get_topography_top_index( 's' ) )
     3581          k_topo_l = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
    35863582#if defined( __parallel )
    35873583          CALL MPI_ALLREDUCE( k_topo_l, k_topo, 1, MPI_INTEGER, MPI_MIN, &
     
    37733769             DO  i = nxl, nxr
    37743770                DO  j = nys, nyn
    3775                    k_topo_l = get_topography_top_index_ji( j, i, 's' )
     3771                   k_topo_l = topo_top_ind(j,i,0)
    37763772                   DO k = k_topo_l+1, nzt+1
    37773773                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo_l)  * d_hours_day
     
    40664062!
    40674063!--             Obtain topography top index (lower bound of RRTMG)
    4068                 k_topo = get_topography_top_index_ji( j, i, 's' )
     4064                k_topo = topo_top_ind(j,i,0)
    40694065
    40704066                IF ( lw_radiation )  THEN
     
    56885684!
    56895685!--          Following expression equals former kk = k - nzb_s_inner(j,i)
    5690              kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
     5686             kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
    56915687             pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
    56925688                 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
     
    57015697                 j = pcbl(iy, ipcgb)
    57025698                 k = pcbl(iz, ipcgb)
    5703                  kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
     5699                 kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
    57045700                 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), &
    57055701                                                   pc_transpiration_rate(kk,j,i), pc_latent_rate(kk,j,i) )
     
    62176213!--    removed later). The following contruct finds the lowest / largest index
    62186214!--    for any upward-facing wall (see bit 12).
    6219        nzubl = MINVAL( get_topography_top_index( 's' ) )
    6220        nzutl = MAXVAL( get_topography_top_index( 's' ) )
     6215       nzubl = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
     6216       nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
    62216217
    62226218       nzubl = MAX( nzubl, nzb )
     
    62356231!
    62366232!--                Find topography top index
    6237                    k_topo = get_topography_top_index_ji( j, i, 's' )
     6233                   k_topo = topo_top_ind(j,i,0)
    62386234
    62396235                   DO k = nzt+1, 0, -1
     
    63586354!
    63596355!--                Find topography top index
    6360                    k_topo = get_topography_top_index_ji( j, i, 's' )
     6356                   k_topo = topo_top_ind(j,i,0)
    63616357
    63626358                   DO k = k_topo + 1, pct(j,i)
     
    67936789        ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
    67946790        nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
    6795         nzterrl = get_topography_top_index( 's' )
     6791        nzterrl = topo_top_ind(nys:nyn,nxl:nxr,0)
    67966792        CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
    67976793                            nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
     
    68036799        DEALLOCATE(nzterrl_l)
    68046800#else
    6805         nzterr = RESHAPE( get_topography_top_index( 's' ), (/(nx+1)*(ny+1)/) )
     6801        nzterr = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
    68066802#endif
    68076803        IF ( plant_canopy )  THEN
     
    68846880            DO i = nxl, nxr
    68856881                DO j = nys, nyn
    6886                     k = get_topography_top_index_ji( j, i, 's' )
     6882                    k = topo_top_ind(j,i,0)
    68876883
    68886884                    sub_lad(k:nz_plant_t, j, i) = lad_s(0:nz_plant_t-k, j, i)
     
    1014310139    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
    1014410140
     10141    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
     10142
    1014510143    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
    1014610144
     
    1054010538    LOGICAL      ::  found       !<
    1054110539
     10540    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
     10541
    1054210542    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
    1054310543
     
    1110711107    CHARACTER (LEN=*) ::  variable   !<
    1110811108
     11109    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
     11110
    1110911111    INTEGER(iwp) ::  av              !<
    1111011112    INTEGER(iwp) ::  i               !<
    1111111113    INTEGER(iwp) ::  j               !<
    11112     INTEGER(iwp) ::  k               !<
    11113     INTEGER(iwp) ::  im              !< loop index for masked variables
    11114     INTEGER(iwp) ::  jm              !< loop index for masked variables
    11115     INTEGER(iwp) ::  kk              !<
     11114    INTEGER(iwp) ::  k               !<
    1111611115    INTEGER(iwp) ::  mid             !< masked output running index
    11117     INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
     11116    INTEGER(iwp) ::  topo_top_index  !< k index of highest horizontal surface
    1111811117
    1111911118    LOGICAL ::  found                !< true if output array was found
     
    1112711126    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
    1112811127
     11128
     11129    found    = .TRUE.
     11130    grid     = 's'
    1112911131    resorted = .FALSE.
    11130     found = .TRUE.
    1113111132
    1113211133    SELECT CASE ( TRIM( variable ) )
     
    1121511216             DO  j = 1, mask_size_l(mid,2)
    1121611217!
    11217 !--             Get k index of the highest terraing surface
    11218                 im = mask_i(mid,i)
    11219                 jm = mask_j(mid,j)
    11220                 ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), &
    11221                               DIM = 1 ) - 1
     11218!--             Get k index of highest horizontal surface
     11219                topo_top_index = topo_top_ind(mask_j(mid,j), &
     11220                                              mask_i(mid,i),   &
     11221                                              0 )
     11222!
     11223!--             Save output array
    1122211224                DO  k = 1, mask_size_l(mid,3)
    11223                    kk = MIN( ktt+mask_k(mid,k), nzt+1 )
    11224 !
    11225 !--                Set value if not in building
    11226                    IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
    11227                       local_pf(i,j,k) = fill_value
    11228                    ELSE
    11229                       local_pf(i,j,k) = to_be_resorted(kk,jm,im)
    11230                    ENDIF
     11225                   local_pf(i,j,k) = to_be_resorted(                         &
     11226                                          MIN( topo_top_index+mask_k(mid,k), &
     11227                                               nzt+1 ),                      &
     11228                                          mask_j(mid,j),                     &
     11229                                          mask_i(mid,i)                     )
    1123111230                ENDDO
    1123211231             ENDDO
Note: See TracChangeset for help on using the changeset viewer.