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

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

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • 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.