Changeset 3129 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jul 16, 2018 7:45:13 AM (6 years ago)
Author:
gronemeier
Message:

merge with branch rans: update of rans mode and data output

Location:
palm/trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/Makefile

  • palm/trunk/SOURCE/boundary_conds.f90

    r2938 r3129  
    2525! -----------------
    2626! $Id$
     27! - Use wall function for e_p and diss_p in case of rans_tke_e
     28! - move limitation of diss_p from tcm_prognostic
     29!
     30! 2938 2018-03-27 15:52:42Z suehring
    2731! Set boundary condition for TKE and TKE dissipation rate in case of nesting
    2832! and if parent model operates in RANS mode but child model in LES mode.
     
    182186    USE arrays_3d,                                                             &
    183187        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
    184                diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s,  &
     188               diss, diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s,  &
    185189               s_p, sa, sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,  &
    186190               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
    187                w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
     191               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init, ddzu
    188192
    189193    USE chemistry_model_mod,                                                   &
     
    196200               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
    197201               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
    198                inflow_s, intermediate_timestep_count,                          &
     202               inflow_s, intermediate_timestep_count, kappa,                   &
    199203               microphysics_morrison, microphysics_seifert, nest_domain,       &
    200204               nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, nudging,&
     
    217221
    218222    USE surface_mod,                                                           &
    219         ONLY :  bc_h
     223        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
     224                surf_usm_h, surf_usm_v
     225
     226    USE turbulence_closure_mod,                                                &
     227        ONLY:  c_0
    220228
    221229    IMPLICIT NONE
     
    344352             ENDDO
    345353          ENDDO
     354       ELSE
     355!
     356!--       Use wall function within constant-flux layer
     357!--       Upward-facing surfaces
     358!--       Default surfaces
     359          DO  m = 1, surf_def_h(0)%ns
     360             i = surf_def_h(0)%i(m)
     361             j = surf_def_h(0)%j(m)
     362             k = surf_def_h(0)%k(m)
     363             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
     364          ENDDO
     365!
     366!--       Natural surfaces
     367          DO  m = 1, surf_lsm_h%ns
     368             i = surf_lsm_h%i(m)
     369             j = surf_lsm_h%j(m)
     370             k = surf_lsm_h%k(m)
     371             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
     372          ENDDO
     373!
     374!--       Urban surfaces
     375          DO  m = 1, surf_usm_h%ns
     376             i = surf_usm_h%i(m)
     377             j = surf_usm_h%j(m)
     378             k = surf_usm_h%k(m)
     379             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
     380          ENDDO
     381!
     382!--       Vertical surfaces
     383          DO  l = 0, 3
     384!
     385!--          Default surfaces
     386             DO  m = 1, surf_def_v(l)%ns
     387                i = surf_def_v(l)%i(m)
     388                j = surf_def_v(l)%j(m)
     389                k = surf_def_v(l)%k(m)
     390                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
     391             ENDDO
     392!
     393!--          Natural surfaces
     394             DO  m = 1, surf_lsm_v(l)%ns
     395                i = surf_lsm_v(l)%i(m)
     396                j = surf_lsm_v(l)%j(m)
     397                k = surf_lsm_v(l)%k(m)
     398                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
     399             ENDDO
     400!
     401!--          Urban surfaces
     402             DO  m = 1, surf_usm_v(l)%ns
     403                i = surf_usm_v(l)%i(m)
     404                j = surf_usm_v(l)%j(m)
     405                k = surf_usm_v(l)%k(m)
     406                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
     407             ENDDO
     408          ENDDO
    346409       ENDIF
    347410
     
    372435!
    373436!-- Boundary conditions for TKE dissipation rate.
    374     IF ( rans_tke_e .AND. .NOT. nest_domain )  THEN
    375        diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
     437    IF ( rans_tke_e )  THEN
     438!
     439!--    Use wall function within constant-flux layer
     440!--    Upward-facing surfaces
     441!--    Default surfaces
     442       DO  m = 1, surf_def_h(0)%ns
     443          i = surf_def_h(0)%i(m)
     444          j = surf_def_h(0)%j(m)
     445          k = surf_def_h(0)%k(m)
     446          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
     447                        / ( kappa * surf_def_h(0)%z_mo(m) )
     448       ENDDO
     449!
     450!--    Natural surfaces
     451       DO  m = 1, surf_lsm_h%ns
     452          i = surf_lsm_h%i(m)
     453          j = surf_lsm_h%j(m)
     454          k = surf_lsm_h%k(m)
     455          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
     456                        / ( kappa * surf_lsm_h%z_mo(m) )
     457       ENDDO
     458!
     459!--    Urban surfaces
     460       DO  m = 1, surf_usm_h%ns
     461          i = surf_usm_h%i(m)
     462          j = surf_usm_h%j(m)
     463          k = surf_usm_h%k(m)
     464          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
     465                        / ( kappa * surf_usm_h%z_mo(m) )
     466       ENDDO
     467!
     468!--    Vertical surfaces
     469       DO  l = 0, 3
     470!
     471!--       Default surfaces
     472          DO  m = 1, surf_def_v(l)%ns
     473             i = surf_def_v(l)%i(m)
     474             j = surf_def_v(l)%j(m)
     475             k = surf_def_v(l)%k(m)
     476             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
     477                           / ( kappa * surf_def_v(l)%z_mo(m) )
     478          ENDDO
     479!
     480!--       Natural surfaces
     481          DO  m = 1, surf_lsm_v(l)%ns
     482             i = surf_lsm_v(l)%i(m)
     483             j = surf_lsm_v(l)%j(m)
     484             k = surf_lsm_v(l)%k(m)
     485             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
     486                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
     487          ENDDO
     488!
     489!--       Urban surfaces
     490          DO  m = 1, surf_usm_v(l)%ns
     491             i = surf_usm_v(l)%i(m)
     492             j = surf_usm_v(l)%j(m)
     493             k = surf_usm_v(l)%k(m)
     494             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
     495                           / ( kappa * surf_usm_v(l)%z_mo(m) )
     496          ENDDO
     497       ENDDO
     498!
     499!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
     500!--    minimum value
     501       DO  i = nxl, nxr
     502          DO  j = nys, nyn
     503             DO  k = nzb, nzt+1
     504                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
     505                                          2.0_wp * diss(k,j,i) ), &
     506                                     0.1_wp * diss(k,j,i),        &
     507                                     0.0001_wp )
     508             ENDDO
     509          ENDDO
     510       ENDDO
     511
     512       IF ( .NOT. nest_domain )  THEN
     513          diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
     514       ENDIF
    376515    ENDIF
    377516
  • palm/trunk/SOURCE/check_parameters.f90

  • palm/trunk/SOURCE/init_3d_model.f90

  • palm/trunk/SOURCE/modules.f90

    r3120 r3129  
    2525! -----------------
    2626! $Id$
     27! add target attribute to km and kh, necessary for output in tcm_data_output_3d
     28!
     29! 3120 2018-07-11 18:30:57Z gronemeier
    2730! +les_dynamic
    2831!
     
    747750    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_v    !< 6th-order advective flux at south face of grid box - v-component
    748751    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_w    !< 6th-order advective flux at south face of grid box - w-component
    749     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kh          !< eddy diffusivity for heat
    750     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  km          !< eddy diffusivity for momentum
     752    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  kh  !< eddy diffusivity for heat
     753    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  km  !< eddy diffusivity for momentum
    751754    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  prr         !< rain rate
    752755    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p_loc       !< local array in multigrid/sor solver containing the pressure which is iteratively advanced in each iteration step
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3121 r3129  
    2525! -----------------
    2626! $Id$
     27! - move limitation of diss to boundary_conds
     28! - move boundary conditions for e and diss to boundary_conds
     29! - consider non-default surfaces in tcm_diffusivities
     30! - use z_mo within surface layer instead of calculating it
     31! - resort output after case select -> reduced code duplication
     32! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
     33!
     34! 3121 2018-07-11 18:46:49Z gronemeier
    2735! - created the function phi_m
    2836! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
     
    117125!> @todo Check for random disturbances
    118126!> @note <Enter notes on the module>
    119 !> @bug  TKE-e closure still crashes due to too small dt
    120127!------------------------------------------------------------------------------!
    121128 MODULE turbulence_closure_mod
     
    202209
    203210    !> @todo remove debug variables
    204     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, &
    205                                                diss_prod2, diss_adve2, diss_diff2, &
    206                                                diss_prod3, diss_adve3, diss_diff3, &
    207                                                dummy1, dummy2, dummy3
     211    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
     212       diss_prod1, diss_adve1, diss_diff1, &
     213       diss_prod2, diss_adve2, diss_diff2, &
     214       diss_prod3, diss_adve3, diss_diff3, &
     215       dummy1, dummy2, dummy3
    208216
    209217
     
    689697    CHARACTER (LEN=*) ::  variable   !< name of variable
    690698
    691     INTEGER(iwp) ::  av   !< flag for (non-)average output
    692     INTEGER(iwp) ::  i    !< loop index
    693     INTEGER(iwp) ::  j    !< loop index
    694     INTEGER(iwp) ::  k    !< loop index
    695     INTEGER(iwp) ::  nzb_do   !< vertical output index (bottom)
    696     INTEGER(iwp) ::  nzt_do   !< vertical output index (top)
     699    INTEGER(iwp) ::  av        !< flag for (non-)average output
     700    INTEGER(iwp) ::  flag_nr   !< number of masking flag
     701    INTEGER(iwp) ::  i         !< loop index
     702    INTEGER(iwp) ::  j         !< loop index
     703    INTEGER(iwp) ::  k         !< loop index
     704    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
     705    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
    697706
    698707    LOGICAL ::  found   !< flag if output variable is found
    699708    LOGICAL ::  two_d   !< flag parameter that indicates 2D variables (horizontal cross sections)
     709    LOGICAL ::  resorted  !< flag if output is already resorted
    700710
    701711    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
     
    704714       !< array to which output data is resorted to
    705715
     716    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
     717   
    706718    found = .TRUE.
     719    resorted = .FALSE.
     720!
     721!-- Set masking flag for topography for not resorted arrays
     722    flag_nr = 0
    707723
    708724    SELECT CASE ( TRIM( variable ) )
    709 
    710725
    711726       CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
    712727          IF ( av == 0 )  THEN
    713              DO  i = nxl, nxr
    714                 DO  j = nys, nyn
    715                    DO k = nzb_do, nzt_do
    716                       local_pf(i,j,k) = diss(k,j,i)
    717                    ENDDO
    718                 ENDDO
    719              ENDDO
     728             to_be_resorted => diss
    720729          ELSE
    721730             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
     
    723732                diss_av = REAL( fill_value, KIND = wp )
    724733             ENDIF
    725              DO  i = nxl, nxr
    726                 DO  j = nys, nyn
    727                    DO k = nzb_do, nzt_do
    728                       local_pf(i,j,k) = diss_av(k,j,i)
    729                    ENDDO
    730                 ENDDO
    731              ENDDO
     734             to_be_resorted => diss_av
    732735          ENDIF
    733 
    734736          IF ( mode == 'xy' ) grid = 'zu'
    735737
    736738       CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
    737739          IF ( av == 0 )  THEN
    738              DO  i = nxl, nxr
    739                 DO  j = nys, nyn
    740                    DO k = nzb_do, nzt_do
    741                       local_pf(i,j,k) = kh(k,j,i)
    742                    ENDDO
    743                 ENDDO
     740             to_be_resorted => kh
     741          ELSE
     742             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
     743                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     744                kh_av = REAL( fill_value, KIND = wp )
     745             ENDIF
     746             to_be_resorted => kh_av
     747          ENDIF
     748          IF ( mode == 'xy' ) grid = 'zu'
     749
     750       CASE ( 'km_xy', 'km_xz', 'km_yz' )
     751          IF ( av == 0 )  THEN
     752             to_be_resorted => km
     753          ELSE
     754             IF ( .NOT. ALLOCATED( km_av ) ) THEN
     755                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     756                km_av = REAL( fill_value, KIND = wp )
     757             ENDIF
     758             to_be_resorted => km_av
     759          ENDIF
     760          IF ( mode == 'xy' ) grid = 'zu'
     761
     762       CASE DEFAULT
     763          found = .FALSE.
     764          grid  = 'none'
     765
     766    END SELECT
     767
     768    IF ( found .AND. .NOT. resorted )  THEN
     769       DO  i = nxl, nxr
     770          DO  j = nys, nyn
     771             DO  k = nzb_do, nzt_do
     772                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
     773                                         REAL( fill_value, KIND = wp ),        &
     774                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
    744775             ENDDO
     776          ENDDO
     777       ENDDO
     778    ENDIF
     779 
     780 END SUBROUTINE tcm_data_output_2d
     781
     782 
     783!------------------------------------------------------------------------------!
     784! Description:
     785! ------------
     786!> Define 3D output variables.
     787!------------------------------------------------------------------------------!
     788 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
     789 
     790
     791    USE averaging,                                                             &
     792        ONLY:  diss_av, kh_av, km_av
     793
     794    IMPLICIT NONE
     795
     796    CHARACTER (LEN=*) ::  variable   !< name of variable
     797
     798    INTEGER(iwp) ::  av        !< flag for (non-)average output
     799    INTEGER(iwp) ::  flag_nr   !< number of masking flag
     800    INTEGER(iwp) ::  i         !< loop index
     801    INTEGER(iwp) ::  j         !< loop index
     802    INTEGER(iwp) ::  k         !< loop index
     803    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
     804    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
     805
     806    LOGICAL ::  found     !< flag if output variable is found
     807    LOGICAL ::  resorted  !< flag if output is already resorted
     808
     809    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
     810
     811    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
     812       !< array to which output data is resorted to
     813
     814    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
     815
     816    found = .TRUE.
     817    resorted = .FALSE.
     818!
     819!-- Set masking flag for topography for not resorted arrays
     820    flag_nr = 0
     821
     822    SELECT CASE ( TRIM( variable ) )
     823
     824       CASE ( 'diss' )
     825          IF ( av == 0 )  THEN
     826             to_be_resorted => diss
    745827          ELSE
    746828             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
     
    748830                diss_av = REAL( fill_value, KIND = wp )
    749831             ENDIF
    750              DO  i = nxl, nxr
    751                 DO  j = nys, nyn
    752                    DO k = nzb_do, nzt_do
    753                       local_pf(i,j,k) = kh_av(k,j,i)
    754                    ENDDO
    755                 ENDDO
    756              ENDDO
    757           ENDIF
    758 
    759           IF ( mode == 'xy' ) grid = 'zu'
    760 
    761        CASE ( 'km_xy', 'km_xz', 'km_yz' )
    762           IF ( av == 0 )  THEN
    763              DO  i = nxl, nxr
    764                 DO  j = nys, nyn
    765                    DO k = nzb_do, nzt_do
    766                       local_pf(i,j,k) = km(k,j,i)
    767                    ENDDO
    768                 ENDDO
    769              ENDDO
    770           ELSE
    771              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
    772                 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    773                 diss_av = REAL( fill_value, KIND = wp )
    774              ENDIF
    775              DO  i = nxl, nxr
    776                 DO  j = nys, nyn
    777                    DO k = nzb_do, nzt_do
    778                       local_pf(i,j,k) = km_av(k,j,i)
    779                    ENDDO
    780                 ENDDO
    781              ENDDO
    782           ENDIF
    783 
    784           IF ( mode == 'xy' ) grid = 'zu'
    785 
    786        CASE DEFAULT
    787           found = .FALSE.
    788           grid  = 'none'
    789 
    790     END SELECT
    791  
    792  END SUBROUTINE tcm_data_output_2d
    793 
    794  
    795 !------------------------------------------------------------------------------!
    796 ! Description:
    797 ! ------------
    798 !> Define 3D output variables.
    799 !------------------------------------------------------------------------------!
    800  SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    801  
    802 
    803     USE averaging,                                                             &
    804         ONLY:  diss_av, kh_av, km_av
    805 
    806     IMPLICIT NONE
    807 
    808     CHARACTER (LEN=*) ::  variable   !< name of variable
    809 
    810     INTEGER(iwp) ::  av     !< flag for (non-)average output
    811     INTEGER(iwp) ::  i      !< loop index
    812     INTEGER(iwp) ::  j      !< loop index
    813     INTEGER(iwp) ::  k      !< loop index
    814     INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
    815     INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
    816 
    817     LOGICAL ::  found   !< flag if output variable is found
    818 
    819     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
    820 
    821     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
    822        !< array to which output data is resorted to
    823 
    824 
    825     found = .TRUE.
    826 
    827 
    828     SELECT CASE ( TRIM( variable ) )
    829 
    830 
    831        CASE ( 'diss' )
    832           IF ( av == 0 )  THEN
    833              DO  i = nxl, nxr
    834                 DO  j = nys, nyn
    835                    DO  k = nzb_do, nzt_do
    836                       local_pf(i,j,k) = diss(k,j,i)
    837                    ENDDO
    838                 ENDDO
    839              ENDDO
    840           ELSE
    841              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
    842                 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    843                 diss_av = REAL( fill_value, KIND = wp )
    844              ENDIF
    845              DO  i = nxl, nxr
    846                 DO  j = nys, nyn
    847                    DO  k = nzb_do, nzt_do
    848                       local_pf(i,j,k) = diss_av(k,j,i)
    849                    ENDDO
    850                 ENDDO
    851              ENDDO
     832             to_be_resorted => diss_av
    852833          ENDIF
    853834
    854835       CASE ( 'kh' )
    855836          IF ( av == 0 )  THEN
    856              DO  i = nxl, nxr
    857                 DO  j = nys, nyn
    858                    DO  k = nzb_do, nzt_do
    859                       local_pf(i,j,k) = kh(k,j,i)
    860                    ENDDO
    861                 ENDDO
    862              ENDDO
     837             to_be_resorted => kh
    863838          ELSE
    864839             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
     
    866841                kh_av = REAL( fill_value, KIND = wp )
    867842             ENDIF
    868              DO  i = nxl, nxr
    869                 DO  j = nys, nyn
    870                    DO  k = nzb_do, nzt_do
    871                       local_pf(i,j,k) = kh_av(k,j,i)
    872                    ENDDO
    873                 ENDDO
    874              ENDDO
     843             to_be_resorted => kh_av
    875844          ENDIF
    876845
    877846       CASE ( 'km' )
    878847          IF ( av == 0 )  THEN
    879              DO  i = nxl, nxr
    880                 DO  j = nys, nyn
    881                    DO  k = nzb_do, nzt_do
    882                       local_pf(i,j,k) = km(k,j,i)
    883                    ENDDO
    884                 ENDDO
    885              ENDDO
     848             to_be_resorted => km
    886849          ELSE
    887850             IF ( .NOT. ALLOCATED( km_av ) ) THEN
     
    889852                km_av = REAL( fill_value, KIND = wp )
    890853             ENDIF
    891              DO  i = nxl, nxr
    892                 DO  j = nys, nyn
    893                    DO  k = nzb_do, nzt_do
    894                       local_pf(i,j,k) = km_av(k,j,i)
    895                    ENDDO
    896                 ENDDO
    897              ENDDO
     854             to_be_resorted => km_av
    898855          ENDIF
    899856
    900857       CASE ( 'dummy3' )                                                        !> @todo remove later
    901858          IF ( av == 0 )  THEN
    902              DO  i = nxl, nxr
    903                 DO  j = nys, nyn
    904                    DO  k = nzb_do, nzt_do
    905                       local_pf(i,j,k) = dummy3(k,j,i)
    906                    ENDDO
    907                 ENDDO
    908              ENDDO
     859             to_be_resorted => dummy3
    909860          ENDIF
    910861
    911862       CASE ( 'diss1' )                                                         !> @todo remove later
    912863          IF ( av == 0 )  THEN
    913              DO  i = nxl, nxr
    914                 DO  j = nys, nyn
    915                    DO  k = nzb_do, nzt_do
    916                       local_pf(i,j,k) = dummy1(k,j,i)
    917                    ENDDO
    918                 ENDDO
    919              ENDDO
     864             to_be_resorted => dummy1
    920865          ENDIF
    921866
    922867       CASE ( 'diss2' )                                                         !> @todo remove later
    923868          IF ( av == 0 )  THEN
    924              DO  i = nxl, nxr
    925                 DO  j = nys, nyn
    926                    DO  k = nzb_do, nzt_do
    927                       local_pf(i,j,k) = dummy2(k,j,i)
    928                    ENDDO
    929                 ENDDO
    930              ENDDO
     869             to_be_resorted => dummy2
    931870          ENDIF
    932871
    933872       CASE ( 'diss_prod1' )                                                    !> @todo remove later
    934873          IF ( av == 0 )  THEN
    935              DO  i = nxl, nxr
    936                 DO  j = nys, nyn
    937                    DO  k = nzb_do, nzt_do
    938                       local_pf(i,j,k) = diss_prod1(k,j,i)
    939                    ENDDO
    940                 ENDDO
    941              ENDDO
     874             to_be_resorted => diss_prod1
    942875          ENDIF
    943876
    944877       CASE ( 'diss_adve1' )                                                    !> @todo remove later
    945878          IF ( av == 0 )  THEN
    946              DO  i = nxl, nxr
    947                 DO  j = nys, nyn
    948                    DO  k = nzb_do, nzt_do
    949                       local_pf(i,j,k) = diss_adve1(k,j,i)
    950                    ENDDO
    951                 ENDDO
    952              ENDDO
     879             to_be_resorted => diss_adve1
    953880          ENDIF
    954881
    955882       CASE ( 'diss_diff1' )                                                    !> @todo remove later
    956883          IF ( av == 0 )  THEN
    957              DO  i = nxl, nxr
    958                 DO  j = nys, nyn
    959                    DO  k = nzb_do, nzt_do
    960                       local_pf(i,j,k) = diss_diff1(k,j,i)
    961                    ENDDO
    962                 ENDDO
    963              ENDDO
     884             to_be_resorted => diss_diff1
    964885          ENDIF
    965886
    966887       CASE ( 'diss_prod2' )                                                    !> @todo remove later
    967888          IF ( av == 0 )  THEN
    968              DO  i = nxl, nxr
    969                 DO  j = nys, nyn
    970                    DO  k = nzb_do, nzt_do
    971                       local_pf(i,j,k) = diss_prod2(k,j,i)
    972                    ENDDO
    973                 ENDDO
    974              ENDDO
     889             to_be_resorted => diss_prod2
    975890          ENDIF
    976891
    977892       CASE ( 'diss_adve2' )                                                    !> @todo remove later
    978893          IF ( av == 0 )  THEN
    979              DO  i = nxl, nxr
    980                 DO  j = nys, nyn
    981                    DO  k = nzb_do, nzt_do
    982                       local_pf(i,j,k) = diss_adve2(k,j,i)
    983                    ENDDO
    984                 ENDDO
    985              ENDDO
     894             to_be_resorted => diss_adve2
    986895          ENDIF
    987896
    988897       CASE ( 'diss_diff2' )                                                    !> @todo remove later
    989898          IF ( av == 0 )  THEN
    990              DO  i = nxl, nxr
    991                 DO  j = nys, nyn
    992                    DO  k = nzb_do, nzt_do
    993                       local_pf(i,j,k) = diss_diff2(k,j,i)
    994                    ENDDO
    995                 ENDDO
    996              ENDDO
     899             to_be_resorted => diss_diff2
    997900          ENDIF
    998901
    999902       CASE ( 'diss_prod3' )                                                    !> @todo remove later
    1000903          IF ( av == 0 )  THEN
    1001              DO  i = nxl, nxr
    1002                 DO  j = nys, nyn
    1003                    DO  k = nzb_do, nzt_do
    1004                       local_pf(i,j,k) = diss_prod3(k,j,i)
    1005                    ENDDO
    1006                 ENDDO
    1007              ENDDO
     904             to_be_resorted => diss_prod3
    1008905          ENDIF
    1009906
    1010907       CASE ( 'diss_adve3' )                                                    !> @todo remove later
    1011908          IF ( av == 0 )  THEN
    1012              DO  i = nxl, nxr
    1013                 DO  j = nys, nyn
    1014                    DO  k = nzb_do, nzt_do
    1015                       local_pf(i,j,k) = diss_adve3(k,j,i)
    1016                    ENDDO
    1017                 ENDDO
    1018              ENDDO
     909             to_be_resorted => diss_adve3
    1019910          ENDIF
    1020911
    1021912       CASE ( 'diss_diff3' )                                                    !> @todo remove later
    1022913          IF ( av == 0 )  THEN
    1023              DO  i = nxl, nxr
    1024                 DO  j = nys, nyn
    1025                    DO  k = nzb_do, nzt_do
    1026                       local_pf(i,j,k) = diss_diff3(k,j,i)
    1027                    ENDDO
    1028                 ENDDO
    1029              ENDDO
     914             to_be_resorted => diss_diff3
    1030915          ENDIF
    1031916         
     
    1034919
    1035920    END SELECT
     921
     922
     923    IF ( found .AND. .NOT. resorted )  THEN
     924       DO  i = nxl, nxr
     925          DO  j = nys, nyn
     926             DO  k = nzb_do, nzt_do
     927                local_pf(i,j,k) = MERGE(                                 &
     928                                   to_be_resorted(k,j,i),                &
     929                                   REAL( fill_value, KIND = wp ),        &
     930                                   BTEST( wall_flags_0(k,j,i), flag_nr ) )
     931             ENDDO
     932          ENDDO
     933       ENDDO
     934       resorted = .TRUE.
     935    ENDIF
    1036936
    1037937 END SUBROUTINE tcm_data_output_3d
     
    11681068
    11691069       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    1170 !
    1171 !--       Transfer initial profiles to the arrays of the 3D model
    1172           DO  i = nxlg, nxrg
    1173              DO  j = nysg, nyng
    1174                 e(:,j,i)  = e1d
    1175                 kh(:,j,i) = kh1d
    1176                 km(:,j,i) = km1d
     1070
     1071          IF ( .NOT. rans_tke_e ) THEN
     1072!
     1073!--          Transfer initial profiles to the arrays of the 3D model
     1074             DO  i = nxlg, nxrg
     1075                DO  j = nysg, nyng
     1076                   e(:,j,i)  = e1d
     1077                   kh(:,j,i) = kh1d
     1078                   km(:,j,i) = km1d
     1079                ENDDO
    11771080             ENDDO
    1178           ENDDO
    1179 
    1180           IF ( constant_diffusion )  THEN
    1181              e = 0.0_wp
    1182           ENDIF
    1183 
    1184           IF ( rans_tke_e )  THEN
    1185              IF ( dissipation_1d == 'prognostic' )  THEN    !> @query Why must this be checked?
    1186                 DO  i = nxlg, nxrg                          !>   Should 'diss' not always
    1187                    DO  j = nysg, nyng                       !>   be prognostic in case rans_tke_e?
    1188                       diss(:,j,i) = diss1d
    1189                    ENDDO
    1190                 ENDDO
    1191              ELSE
     1081
     1082             IF ( constant_diffusion )  THEN
     1083                e = 0.0_wp
     1084             ENDIF
     1085
     1086          ELSE
     1087!
     1088!--          In case of TKE-e closure in RANS mode, do not use e, diss, and km
     1089!--          profiles from 1D model. Instead, initialize with constant profiles
     1090             IF ( constant_diffusion )  THEN
     1091                km = km_constant
     1092                kh = km / prandtl_number
     1093                e  = 0.0_wp
     1094             ELSEIF ( e_init > 0.0_wp )  THEN
    11921095                DO  i = nxlg, nxrg
    11931096                   DO  j = nysg, nyng
    11941097                      DO  k = nzb+1, nzt
    1195                          diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k)
     1098                         km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
    11961099                      ENDDO
    11971100                   ENDDO
    11981101                ENDDO
     1102                km(nzb,:,:)   = km(nzb+1,:,:)
     1103                km(nzt+1,:,:) = km(nzt,:,:)
     1104                kh = km / prandtl_number
     1105                e  = e_init
     1106             ELSE
     1107                IF ( .NOT. ocean )  THEN
     1108                   kh   = 0.01_wp   ! there must exist an initial diffusion, because
     1109                   km   = 0.01_wp   ! otherwise no TKE would be produced by the
     1110                                     ! production terms, as long as not yet
     1111                                     ! e = (u*/cm)**2 at k=nzb+1
     1112                ELSE
     1113                   kh   = 0.00001_wp
     1114                   km   = 0.00001_wp
     1115                ENDIF
     1116                e    = 0.0_wp
    11991117             ENDIF
     1118
     1119             DO  i = nxlg, nxrg
     1120                DO  j = nysg, nyng
     1121                   DO  k = nzb+1, nzt
     1122                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
     1123                   ENDDO
     1124                ENDDO
     1125             ENDDO
     1126             diss(nzb,:,:) = diss(nzb+1,:,:)
     1127             diss(nzt+1,:,:) = diss(nzt,:,:)
     1128
    12001129          ENDIF
    12011130
     
    13811310
    13821311    USE control_parameters,                                                    &
    1383         ONLY:  bc_lr_cyc, bc_ns_cyc, f, kappa, message_string,                 &
    1384                wall_adjustment_factor
     1312        ONLY:  bc_lr_cyc, bc_ns_cyc, f, message_string, wall_adjustment_factor
    13851313
    13861314    USE grid_variables,                                                        &
     
    19731901
    19741902    USE surface_mod,                                                           &
    1975         ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h
     1903        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
    19761904
    19771905    IMPLICIT NONE
     
    25562484
    25572485!
    2558 !--    Use special boundary condition in case of TKE-e closure
    2559        IF ( rans_tke_e )  THEN
    2560           surf_s = surf_def_h(0)%start_index(j,i)
    2561           surf_e = surf_def_h(0)%end_index(j,i)
    2562           DO  m = surf_s, surf_e
    2563              k = surf_def_h(0)%k(m)
    2564              e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
    2565           ENDDO
    2566 
    2567           DO  l = 0, 3
    2568              surf_s = surf_def_v(l)%start_index(j,i)
    2569              surf_e = surf_def_v(l)%end_index(j,i)
    2570              DO  m = surf_s, surf_e
    2571                 k = surf_def_v(l)%k(m)
    2572                 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
    2573              ENDDO
    2574           ENDDO
    2575        ENDIF
    2576 
    2577 !
    25782486!--    Calculate tendencies for the next Runge-Kutta step
    25792487       IF ( timestep_scheme(1:5) == 'runge' )  THEN
     
    26752583
    26762584!
    2677 !--    Use special boundary condition in case of TKE-e closure
    2678        surf_s = surf_def_h(0)%start_index(j,i)
    2679        surf_e = surf_def_h(0)%end_index(j,i)
    2680        DO  m = surf_s, surf_e
    2681           k = surf_def_h(0)%k(m)
    2682           diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / kappa * ddzu(k)
    2683        ENDDO
    2684 
    2685        DO  l = 0, 1
    2686           surf_s = surf_def_v(l)%start_index(j,i)
    2687           surf_e = surf_def_v(l)%end_index(j,i)
    2688           DO  m = surf_s, surf_e
    2689              k = surf_def_v(l)%k(m)
    2690              diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy )
    2691           ENDDO
    2692        ENDDO
    2693 
    2694        DO  l = 2, 3
    2695           surf_s = surf_def_v(l)%start_index(j,i)
    2696           surf_e = surf_def_v(l)%end_index(j,i)
    2697           DO  m = surf_s, surf_e
    2698              k = surf_def_v(l)%k(m)
    2699              diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx )
    2700           ENDDO
    2701        ENDDO
    2702 !
    27032585!--    Calculate tendencies for the next Runge-Kutta step
    27042586       IF ( timestep_scheme(1:5) == 'runge' )  THEN
     
    27152597          ENDIF
    27162598       ENDIF
    2717 !
    2718 !--    Limit change of diss to be between -90% and +100%. Also, set an absolute
    2719 !--    minimum value
    2720        DO  k = nzb, nzt+1
    2721           diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i),         &
    2722                                     0.1_wp * diss(k,j,i),  &
    2723                                     0.0001_wp ),           &
    2724                                2.0_wp * diss(k,j,i) )
    2725        ENDDO
    27262599
    27272600       IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = diss_p(:,j,i)   !> @todo remove later
     
    47774650!> Computation of the turbulent diffusion coefficients for momentum and heat
    47784651!> according to Prandtl-Kolmogorov.
    4779 !> @todo consider non-default surfaces
    47804652!------------------------------------------------------------------------------!
    47814653 SUBROUTINE tcm_diffusivities_default( var, var_reference )
     
    47924664
    47934665    USE surface_mod,                                                           &
    4794         ONLY :  bc_h, surf_def_h, surf_def_v
     4666        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
     4667                surf_usm_h, surf_usm_v
    47954668
    47964669    IMPLICIT NONE
     
    49564829!
    49574830!--    Up- and downward facing surfaces
     4831!--    Default surfaces
    49584832       DO  n = 0, 1
    49594833          DO  m = 1, surf_def_h(n)%ns
     
    49614835             j = surf_def_h(n)%j(m)
    49624836             k = surf_def_h(n)%k(m)
    4963              km(k,j,i) = kappa * surf_def_h(n)%us(m) * dzu(k)
     4837             km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m)
    49644838             kh(k,j,i) = 1.35_wp * km(k,j,i)
    49654839          ENDDO
    49664840       ENDDO
    49674841!
    4968 !--    North- and southward facing surfaces
    4969        DO  n = 0, 1
     4842!--    Upward facing surfaces
     4843!--    Natural surfaces
     4844       DO  m = 1, surf_lsm_h%ns
     4845          i = surf_lsm_h%i(m)
     4846          j = surf_lsm_h%j(m)
     4847          k = surf_lsm_h%k(m)
     4848          km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m)
     4849          kh(k,j,i) = 1.35_wp * km(k,j,i)
     4850       ENDDO
     4851!
     4852!--    Urban surfaces
     4853       DO  m = 1, surf_usm_h%ns
     4854          i = surf_usm_h%i(m)
     4855          j = surf_usm_h%j(m)
     4856          k = surf_usm_h%k(m)
     4857          km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m)
     4858          kh(k,j,i) = 1.35_wp * km(k,j,i)
     4859       ENDDO
     4860
     4861!
     4862!--    North-, south-, west and eastward facing surfaces
     4863       DO  n = 0, 3
     4864!
     4865!--       Default surfaces
    49704866          DO  m = 1, surf_def_v(n)%ns
    49714867             i = surf_def_v(n)%i(m)
    49724868             j = surf_def_v(n)%j(m)
    49734869             k = surf_def_v(n)%k(m)
    4974              km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy
     4870             km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
    49754871             kh(k,j,i) = 1.35_wp * km(k,j,i)
    49764872          ENDDO
    4977        ENDDO
    4978 !
    4979 !--    West- and eastward facing surfaces
    4980        DO  n = 2, 3
    4981           DO  m = 1, surf_def_v(n)%ns
    4982              i = surf_def_v(n)%i(m)
    4983              j = surf_def_v(n)%j(m)
    4984              k = surf_def_v(n)%k(m)
    4985              km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx
     4873!
     4874!--       Natural surfaces
     4875          DO  m = 1, surf_lsm_v(n)%ns
     4876             i = surf_lsm_v(n)%i(m)
     4877             j = surf_lsm_v(n)%j(m)
     4878             k = surf_lsm_v(n)%k(m)
     4879             km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
     4880             kh(k,j,i) = 1.35_wp * km(k,j,i)
     4881          ENDDO
     4882!
     4883!--       Urban surfaces
     4884          DO  m = 1, surf_usm_v(n)%ns
     4885             i = surf_usm_v(n)%i(m)
     4886             j = surf_usm_v(n)%j(m)
     4887             k = surf_usm_v(n)%k(m)
     4888             km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
    49864889             kh(k,j,i) = 1.35_wp * km(k,j,i)
    49874890          ENDDO
Note: See TracChangeset for help on using the changeset viewer.