Changeset 2968 for palm/trunk/SOURCE


Ignore:
Timestamp:
Apr 13, 2018 11:52:24 AM (7 years ago)
Author:
suehring
Message:

Bugfixes in initalization of land-surface model as well as timeseries data output in case of elevated model surfaces.

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r2817 r2968  
    2525! -----------------
    2626! $Id$
     27! Bugfix in output of timeseries quantities in case of elevated model surfaces.
     28!
     29! 2817 2018-02-19 16:32:21Z knoop
    2730! Preliminary gust module interface implemented
    2831!
     
    292295    USE indices,                                                               &
    293296        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
    294                 ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
     297                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, topo_min_level,     &
     298                wall_flags_0
    295299       
    296300    USE kinds
     
    21022106
    21032107!
    2104 !--    Collect the time series quantities
    2105        ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
    2106        ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
     2108!--    Collect the time series quantities. Please note, timeseries quantities
     2109!--    which are collected from horizontally averaged profiles, e.g. wpt
     2110!--    or pt(zp), are treated specially. In case of elevated model surfaces,
     2111!--    index nzb+1 might be within topography and data will be zero. Therefore,
     2112!--    take value for the first atmosphere index, which is topo_min_level+1.
     2113       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
     2114       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
    21072115       ts_value(3,sr) = dt_3d
    2108        ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
    2109        ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
     2116       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)          ! u*
     2117       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)        ! th*
    21102118       ts_value(6,sr) = u_max
    21112119       ts_value(7,sr) = v_max
    21122120       ts_value(8,sr) = w_max
    2113        ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
    2114        ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
    2115        ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
    2116        ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
    2117        ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
    2118        ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
    2119        ts_value(15,sr) = hom(nzb+1,1,16,sr)        ! w'pt'   at k=1
    2120        ts_value(16,sr) = hom(nzb+1,1,18,sr)        ! wpt     at k=1
    2121        ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)   ! pt(0)
    2122        ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
    2123        ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
    2124        ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
    2125        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
     2121       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)       ! new divergence
     2122       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)       ! old Divergence
     2123       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)       ! z_i(1)
     2124       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)       ! z_i(2)
     2125       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)       ! w*
     2126       ts_value(14,sr) = hom(nzb,1,16,sr)              ! w'pt'   at k=0
     2127       ts_value(15,sr) = hom(topo_min_level+1,1,16,sr) ! w'pt'   at k=1
     2128       ts_value(16,sr) = hom(topo_min_level+1,1,18,sr) ! wpt     at k=1
     2129       ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)      ! pt(0)
     2130       ts_value(18,sr) = hom(topo_min_level+1,1,4,sr)  ! pt(zp)
     2131       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)       ! u'w'    at k=0
     2132       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)       ! v'w'    at k=0
     2133       ts_value(21,sr) = hom(nzb,1,48,sr)              ! w"q"    at k=0
    21262134
    21272135       IF ( .NOT. neutral )  THEN
  • palm/trunk/SOURCE/init_grid.f90

    r2955 r2968  
    2525! -----------------
    2626! $Id$
     27! Bugfix in initialization in case of elevated model surface. Introduce
     28! index for minimum topography-top.
     29!
     30! 2955 2018-04-09 15:14:01Z suehring
    2731! Improve topography filter routine and add ghost-point exchange for building
    2832! ID and building type.
     
    332336               nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,                 &
    333337               nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner,             &
    334                nzb_w_outer, nzt
     338               nzb_w_outer, nzt, topo_min_level
    335339   
    336340    USE kinds
     
    348352    IMPLICIT NONE
    349353
    350     INTEGER(iwp) ::  i             !< index variable along x
    351     INTEGER(iwp) ::  j             !< index variable along y
    352     INTEGER(iwp) ::  k             !< index variable along z
    353     INTEGER(iwp) ::  k_top         !< topography top index on local PE
    354     INTEGER(iwp) ::  l             !< loop variable
    355     INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
    356     INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
     354    INTEGER(iwp) ::  i                !< index variable along x
     355    INTEGER(iwp) ::  j                !< index variable along y
     356    INTEGER(iwp) ::  k                !< index variable along z
     357    INTEGER(iwp) ::  k_top            !< topography top index on local PE
     358    INTEGER(iwp) ::  l                !< loop variable
     359    INTEGER(iwp) ::  nzb_local_max    !< vertical grid index of maximum topography height
     360    INTEGER(iwp) ::  nzb_local_min    !< vertical grid index of minimum topography height
    357361                                     
    358362    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
     
    553557!   
    554558!-- Finally, if topography extents up to the model top, limit nzb_max to nzt.
    555     nzb_max = MIN( nzb_max, nzt )
    556 
     559    nzb_max = MIN( nzb_max, nzt )
     560!
     561!-- Determine minimum index of topography. Usually, this will be nzb. In case
     562!-- there is elevated topography, however, the lowest topography will be higher.
     563!-- This index is e.g. used to calculate mean first-grid point atmosphere
     564!-- temperature, surface pressure and density, etc. .
     565    topo_min_level   = 0
     566#if defined( __parallel )
     567    CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),             &
     568                        topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
     569#else
     570    topo_min_level = MINVAL( get_topography_top_index( 's' ) )
     571#endif
    557572!
    558573!-- Initialize boundary conditions via surface type
    559574    CALL init_bc
    560 
    561575!
    562576!-- Allocate and set topography height arrays required for data output
     
    621635    nzb_local(nys:nyn,nxl:nxr) = get_topography_top_index( 's' )
    622636    CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    623 
     637!
     638!-- Check topography for consistency with model domain. Therefore, use
     639!-- maximum and minium topography-top indices. Note, minimum topography top
     640!-- index is already calculated. 
    624641    IF ( TRIM( topography ) /= 'flat' )  THEN
    625642#if defined( __parallel )
    626643       CALL MPI_ALLREDUCE( MAXVAL( get_topography_top_index( 's' ) ),          &
    627                            nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
    628        CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),          &
    629                            nzb_local_min, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )                   
     644                           nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )             
    630645#else
    631646       nzb_local_max = MAXVAL( get_topography_top_index( 's' ) )
    632        nzb_local_min = MINVAL( get_topography_top_index( 's' ) )
    633647#endif
     648       nzb_local_min = topo_min_level
    634649!
    635650!--    Consistency checks
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2963 r2968  
    2525! -----------------
    2626! $Id$
     27! Bugfix in initialization in case of elevated model surface
     28!
     29! 2963 2018-04-12 14:47:44Z suehring
    2730! - In initialization of surface types, consider the case that surface_fractions
    2831!   is not given in static input file.
     
    23742377
    23752378       USE indices,                                                            &
    2376            ONLY:  nx, ny
     2379           ONLY:  nx, ny, topo_min_level
    23772380
    23782381       USE pmc_interface,                                                      &
     
    24082411       IF (  .NOT.  cloud_physics )  THEN
    24092412          CALL calc_mean_profile( pt, 4 )
    2410           rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn )
     2413          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(topo_min_level+1,1,4,0) * exn )
    24112414       ENDIF
    24122415
     
    24242427       tt_soil_h_m%var_2d    = 0.0_wp
    24252428       tm_soil_h_m%var_2d    = 0.0_wp
    2426        tm_liq_h_m%var_1d  = 0.0_wp
    2427        surf_lsm_h%c_liq = 0.0_wp
     2429       tm_liq_h_m%var_1d     = 0.0_wp
     2430       surf_lsm_h%c_liq      = 0.0_wp
    24282431
    24292432       surf_lsm_h%ghf = 0.0_wp
     
    24432446          tt_soil_v_m(l)%var_2d    = 0.0_wp
    24442447          tm_soil_v_m(l)%var_2d    = 0.0_wp
    2445           tm_liq_v_m(l)%var_1d  = 0.0_wp
    2446           surf_lsm_v(l)%c_liq = 0.0_wp
     2448          tm_liq_v_m(l)%var_1d     = 0.0_wp
     2449          surf_lsm_v(l)%c_liq      = 0.0_wp
    24472450
    24482451          surf_lsm_v(l)%ghf = 0.0_wp
  • palm/trunk/SOURCE/modules.f90

    r2964 r2968  
    2525! -----------------
    2626! $Id$
     27! +topo_min_level
     28!
     29! 2964 2018-04-12 16:04:03Z raasch
    2730! *_time_count variables are all initialized with zero now
    2831!
     
    17191722    USE kinds
    17201723
    1721     INTEGER(iwp) ::  nbgp = 3     !< number of boundary ghost points
    1722     INTEGER(iwp) ::  ngp_sums     !< number of vertical profile grid points time number of output profiles - used for allreduce statements in MPI calls
    1723     INTEGER(iwp) ::  ngp_sums_ls  !< number of vertical profile grid points time number of large-scale forcing profiles - used for allreduce statements in MPI calls
    1724     INTEGER(iwp) ::  nnx          !< number of subdomain grid points in x-direction
    1725     INTEGER(iwp) ::  nx = 0       !< nx+1 = total number of grid points in x-direction
    1726     INTEGER(iwp) ::  nx_a         !< in coupled atmosphere-ocean runs: total number of grid points along x (atmosphere)
    1727     INTEGER(iwp) ::  nx_o         !< in coupled atmosphere-ocean runs: total number of grid points along x (ocean)
    1728     INTEGER(iwp) ::  nxl          !< left-most grid index of subdomain (excluding ghost points)
    1729     INTEGER(iwp) ::  nxlg         !< left-most grid index of subdomain (including ghost points)
    1730     INTEGER(iwp) ::  nxlu         !< =nxl+1 (at left domain boundary with inflow from left), else =nxl (used for u-velocity component)
    1731     INTEGER(iwp) ::  nxr          !< right-most grid index of subdomain (excluding ghost points)
    1732     INTEGER(iwp) ::  nxrg         !< right-most grid index of subdomain (including ghost points)
    1733     INTEGER(iwp) ::  nx_on_file   !< nx of previous run in job chain
    1734     INTEGER(iwp) ::  nny          !< number of subdomain grid points in y-direction
    1735     INTEGER(iwp) ::  ny = 0       !< ny+1 = total number of grid points in y-direction
    1736     INTEGER(iwp) ::  ny_a         !< in coupled atmosphere-ocean runs: total number of grid points along y (atmosphere)
    1737     INTEGER(iwp) ::  ny_o         !< in coupled atmosphere-ocean runs: total number of grid points along y (ocean)
    1738     INTEGER(iwp) ::  nyn          !< north-most grid index of subdomain (excluding ghost points)
    1739     INTEGER(iwp) ::  nyng         !< north-most grid index of subdomain (including ghost points)
    1740     INTEGER(iwp) ::  nys          !< south-most grid index of subdomain (excluding ghost points)
    1741     INTEGER(iwp) ::  nysg         !< south-most grid index of subdomain (including ghost points)
    1742     INTEGER(iwp) ::  nysv         !< =nys+1 (at south domain boundary with inflow from south), else =nys (used for v-velocity component)
    1743     INTEGER(iwp) ::  ny_on_file   !< ny of previous run in job chain
    1744     INTEGER(iwp) ::  nnz          !< number of subdomain grid points in z-direction
    1745     INTEGER(iwp) ::  nz = 0       !< total number of grid points in z-direction
    1746     INTEGER(iwp) ::  nzb          !< bottom grid index of computational domain
    1747     INTEGER(iwp) ::  nzb_diff     !< will be removed
    1748     INTEGER(iwp) ::  nzb_max      !< vertical index of topography top
    1749     INTEGER(iwp) ::  nzt          !< nzt+1 = top grid index of computational domain
     1724    INTEGER(iwp) ::  nbgp = 3       !< number of boundary ghost points
     1725    INTEGER(iwp) ::  ngp_sums       !< number of vertical profile grid points time number of output profiles - used for allreduce statements in MPI calls
     1726    INTEGER(iwp) ::  ngp_sums_ls    !< number of vertical profile grid points time number of large-scale forcing profiles - used for allreduce statements in MPI calls
     1727    INTEGER(iwp) ::  nnx            !< number of subdomain grid points in x-direction
     1728    INTEGER(iwp) ::  nx = 0         !< nx+1 = total number of grid points in x-direction
     1729    INTEGER(iwp) ::  nx_a           !< in coupled atmosphere-ocean runs: total number of grid points along x (atmosphere)
     1730    INTEGER(iwp) ::  nx_o           !< in coupled atmosphere-ocean runs: total number of grid points along x (ocean)
     1731    INTEGER(iwp) ::  nxl            !< left-most grid index of subdomain (excluding ghost points)
     1732    INTEGER(iwp) ::  nxlg           !< left-most grid index of subdomain (including ghost points)
     1733    INTEGER(iwp) ::  nxlu           !< =nxl+1 (at left domain boundary with inflow from left), else =nxl (used for u-velocity component)
     1734    INTEGER(iwp) ::  nxr            !< right-most grid index of subdomain (excluding ghost points)
     1735    INTEGER(iwp) ::  nxrg           !< right-most grid index of subdomain (including ghost points)
     1736    INTEGER(iwp) ::  nx_on_file     !< nx of previous run in job chain
     1737    INTEGER(iwp) ::  nny            !< number of subdomain grid points in y-direction
     1738    INTEGER(iwp) ::  ny = 0         !< ny+1 = total number of grid points in y-direction
     1739    INTEGER(iwp) ::  ny_a           !< in coupled atmosphere-ocean runs: total number of grid points along y (atmosphere)
     1740    INTEGER(iwp) ::  ny_o           !< in coupled atmosphere-ocean runs: total number of grid points along y (ocean)
     1741    INTEGER(iwp) ::  nyn            !< north-most grid index of subdomain (excluding ghost points)
     1742    INTEGER(iwp) ::  nyng           !< north-most grid index of subdomain (including ghost points)
     1743    INTEGER(iwp) ::  nys            !< south-most grid index of subdomain (excluding ghost points)
     1744    INTEGER(iwp) ::  nysg           !< south-most grid index of subdomain (including ghost points)
     1745    INTEGER(iwp) ::  nysv           !< =nys+1 (at south domain boundary with inflow from south), else =nys (used for v-velocity component)
     1746    INTEGER(iwp) ::  ny_on_file     !< ny of previous run in job chain
     1747    INTEGER(iwp) ::  nnz            !< number of subdomain grid points in z-direction
     1748    INTEGER(iwp) ::  nz = 0         !< total number of grid points in z-direction
     1749    INTEGER(iwp) ::  nzb            !< bottom grid index of computational domain
     1750    INTEGER(iwp) ::  nzb_diff       !< will be removed
     1751    INTEGER(iwp) ::  nzb_max        !< vertical index of topography top
     1752    INTEGER(iwp) ::  nzt            !< nzt+1 = top grid index of computational domain
     1753    INTEGER(iwp) ::  topo_min_level !< minimum topography-top index (usually equal to nzb)
    17501754
    17511755    INTEGER(idp), DIMENSION(:), ALLOCATABLE ::  ngp_3d        !< number of grid points of the total domain
Note: See TracChangeset for help on using the changeset viewer.