!> @file turbulence_closure_mod.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 2017-2018 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: turbulence_closure_mod.f90 2764 2018-01-22 09:25:36Z thiele $ ! Bugfix: remove duplicate SAVE statements ! ! 2746 2018-01-15 12:06:04Z suehring ! Move flag plant canopy to modules ! ! 2718 2018-01-02 08:49:38Z maronga ! Corrected "Former revisions" section ! ! 2701 2017-12-15 15:40:50Z suehring ! Changes from last commit documented ! ! 2698 2017-12-14 18:46:24Z suehring ! Bugfix in get_topography_top_index ! ! 2696 2017-12-14 17:12:51Z kanani ! Initial revision ! ! ! ! ! Authors: ! -------- ! @author Tobias Gronemeier ! ! ! Description: ! ------------ !> This module contains the available turbulence closures for PALM. !> !> !> @todo test initialization for all possibilities !> add OpenMP directives whereever possible !> remove debug output variables (dummy1, dummy2, dummy3) !> @note !> @bug TKE-e closure still crashes due to too small dt !------------------------------------------------------------------------------! MODULE turbulence_closure_mod #if defined( __nopointer ) USE arrays_3d, & ONLY: diss, diss_p, dzu, e, e_p, kh, km, l_grid, & mean_inflow_profiles, prho, pt, tdiss_m, te_m, tend, u, v, vpt, w #else USE arrays_3d, & ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, & e_p, kh, km, l_grid, mean_inflow_profiles, prho, pt, tdiss_m, & te_m, tend, u, v, vpt, w #endif USE control_parameters, & ONLY: constant_diffusion, dt_3d, e_init, humidity, inflow_l, & initializing_actions, intermediate_timestep_count, & intermediate_timestep_count_max, kappa, km_constant, les_mw, & ocean, plant_canopy, prandtl_number, prho_reference, & pt_reference, rans_mode, rans_tke_e, rans_tke_l, simulated_time,& timestep_scheme, turbulence_closure, turbulent_inflow, & use_upstream_for_tke, vpt_reference, ws_scheme_sca USE advec_ws, & ONLY: advec_s_ws USE advec_s_bc_mod, & ONLY: advec_s_bc USE advec_s_pw_mod, & ONLY: advec_s_pw USE advec_s_up_mod, & ONLY: advec_s_up USE cpulog, & ONLY: cpu_log, log_point, log_point_s USE indices, & ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, & nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, & wall_flags_0 USE kinds USE pegrid USE plant_canopy_model_mod, & ONLY: pcm_tendency USE statistics, & ONLY: hom, hom_sum, statistic_regions USE user_actions_mod, & ONLY: user_actions IMPLICIT NONE REAL(wp) :: c_1 = 1.44_wp !< model constant for RANS mode REAL(wp) :: c_2 = 1.92_wp !< model constant for RANS mode REAL(wp) :: c_3 = 1.44_wp !< model constant for RANS mode REAL(wp) :: c_h = 0.0015_wp !< model constant for RANS mode REAL(wp) :: c_m !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES) REAL(wp) :: c_mu = 0.09_wp !< model constant for RANS mode REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length REAL(wp) :: sig_e = 1.0_wp !< factor to calculate Ke from Km REAL(wp) :: sig_diss = 1.3_wp !< factor to calculate K_diss from Km INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy1 !< debug output variable REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy2 !< debug output variable REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy3 !< debug output variable PUBLIC c_m, dummy1, dummy2, dummy3 ! !-- PALM interfaces: !-- Input parameter checks to be done in check_parameters INTERFACE tcm_check_parameters MODULE PROCEDURE tcm_check_parameters END INTERFACE tcm_check_parameters ! !-- Data output checks for 2D/3D data to be done in check_parameters INTERFACE tcm_check_data_output MODULE PROCEDURE tcm_check_data_output END INTERFACE tcm_check_data_output ! !-- Definition of data output quantities INTERFACE tcm_define_netcdf_grid MODULE PROCEDURE tcm_define_netcdf_grid END INTERFACE tcm_define_netcdf_grid ! !-- Averaging of 3D data for output INTERFACE tcm_3d_data_averaging MODULE PROCEDURE tcm_3d_data_averaging END INTERFACE tcm_3d_data_averaging ! !-- Data output of 2D quantities INTERFACE tcm_data_output_2d MODULE PROCEDURE tcm_data_output_2d END INTERFACE tcm_data_output_2d ! !-- Data output of 3D data INTERFACE tcm_data_output_3d MODULE PROCEDURE tcm_data_output_3d END INTERFACE tcm_data_output_3d ! !-- Initialization actions INTERFACE tcm_init MODULE PROCEDURE tcm_init END INTERFACE tcm_init ! !-- Initialization of arrays INTERFACE tcm_init_arrays MODULE PROCEDURE tcm_init_arrays END INTERFACE tcm_init_arrays ! !-- Initialization of TKE production term INTERFACE production_e_init MODULE PROCEDURE production_e_init END INTERFACE production_e_init ! !-- Prognostic equations for TKE and TKE dissipation rate INTERFACE tcm_prognostic MODULE PROCEDURE tcm_prognostic MODULE PROCEDURE tcm_prognostic_ij END INTERFACE tcm_prognostic ! !-- Production term for TKE INTERFACE production_e MODULE PROCEDURE production_e MODULE PROCEDURE production_e_ij END INTERFACE production_e ! !-- Diffusion term for TKE INTERFACE diffusion_e MODULE PROCEDURE diffusion_e MODULE PROCEDURE diffusion_e_ij END INTERFACE diffusion_e ! !-- Diffusion term for TKE dissipation rate INTERFACE diffusion_diss MODULE PROCEDURE diffusion_diss MODULE PROCEDURE diffusion_diss_ij END INTERFACE diffusion_diss ! !-- Mixing length for LES case INTERFACE mixing_length_les MODULE PROCEDURE mixing_length_les END INTERFACE mixing_length_les ! !-- Mixing length for RANS case INTERFACE mixing_length_rans MODULE PROCEDURE mixing_length_rans END INTERFACE mixing_length_rans ! !-- Calculate diffusivities INTERFACE tcm_diffusivities MODULE PROCEDURE tcm_diffusivities END INTERFACE tcm_diffusivities ! !-- Swapping of time levels (required for prognostic variables) INTERFACE tcm_swap_timelevel MODULE PROCEDURE tcm_swap_timelevel END INTERFACE tcm_swap_timelevel SAVE PRIVATE ! !-- Add INTERFACES that must be available to other modules (alphabetical order) PUBLIC production_e_init, tcm_3d_data_averaging, tcm_check_data_output, & tcm_check_parameters, tcm_data_output_2d, tcm_data_output_3d, & tcm_define_netcdf_grid, tcm_diffusivities, tcm_init, & tcm_init_arrays, tcm_prognostic, tcm_swap_timelevel CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check parameters routine for turbulence closure module. !------------------------------------------------------------------------------! SUBROUTINE tcm_check_parameters USE control_parameters, & ONLY: message_string, neutral, turbulent_inflow, turbulent_outflow IMPLICIT NONE ! !-- Define which turbulence closure is going to be used IF ( rans_mode ) THEN c_m = 0.4_wp !according to Detering and Etling (1985) SELECT CASE ( TRIM( turbulence_closure ) ) CASE ( 'TKE-l' ) rans_tke_l = .TRUE. CASE ( 'TKE-e' ) rans_tke_e = .TRUE. IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) & == 0 ) THEN message_string = 'Initializing without 1D model while ' // & 'using TKE-e closure&' // & 'is not possible at the moment!' CALL message( 'tcm_check_parameters', 'TG0005', 1, 2, 0, 6, 0 ) ENDIF CASE DEFAULT message_string = 'Unknown turbulence closure: ' // & TRIM( turbulence_closure ) CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 ) END SELECT message_string = 'RANS mode is still in development! ' // & '&Not all features of PALM are yet compatible '// & 'with RANS mode. &Use at own risk!' CALL message( 'tcm_check_parameters', 'TG0003', 0, 1, 0, 6, 0 ) ELSE c_m = 0.1_wp !according to Lilly (1967) and Deardorff (1980) SELECT CASE ( TRIM( turbulence_closure ) ) CASE ( 'Moeng_Wyngaard' ) les_mw = .TRUE. CASE DEFAULT message_string = 'Unknown turbulence closure: ' // & TRIM( turbulence_closure ) CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 ) END SELECT ENDIF IF ( rans_tke_e ) THEN IF ( turbulent_inflow .OR. turbulent_outflow ) THEN message_string = 'turbulent inflow/outflow is not yet '// & 'implemented for TKE-e closure' CALL message( 'tcm_check_parameters', 'TG0002', 1, 2, 0, 6, 0 ) ENDIF ENDIF END SUBROUTINE tcm_check_parameters !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check data output. !------------------------------------------------------------------------------! SUBROUTINE tcm_check_data_output( var, unit, i, ilen, k ) USE control_parameters, & ONLY: data_output, message_string IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: var !< INTEGER(iwp) :: i !< INTEGER(iwp) :: ilen !< INTEGER(iwp) :: k !< SELECT CASE ( TRIM( var ) ) CASE ( 'diss' ) IF ( .NOT. rans_tke_e ) THEN message_string = 'output of "' // TRIM( var ) // '" requi' // & 'res TKE-e closure for RANS mode.' CALL message( 'tcm_check_data_output', 'TG0101', 1, 2, 0, 6, 0 ) ENDIF unit = 'm2/s3' CASE ( 'dummy2', 'dummy3', 'dummy1' ) unit = '?' CASE ( 'kh', 'km' ) unit = 'm2/s' CASE DEFAULT unit = 'illegal' END SELECT END SUBROUTINE tcm_check_data_output !------------------------------------------------------------------------------! ! Description: ! ------------ !> Define appropriate grid for netcdf variables. !> It is called out from subroutine netcdf. !------------------------------------------------------------------------------! SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) IMPLICIT NONE CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< CHARACTER (LEN=*), INTENT(IN) :: var !< LOGICAL, INTENT(OUT) :: found !< found = .TRUE. ! !-- Check for the grid SELECT CASE ( TRIM( var ) ) CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' ) grid_x = 'x' grid_y = 'y' grid_z = 'zu' CASE ( 'dummy2', 'dummy3', 'dummy1' ) !### remove later grid_x = 'x' grid_y = 'y' grid_z = 'zu' CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' ) grid_x = 'x' grid_y = 'y' grid_z = 'zu' CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' ) grid_x = 'x' grid_y = 'y' grid_z = 'zu' CASE DEFAULT found = .FALSE. grid_x = 'none' grid_y = 'none' grid_z = 'none' END SELECT END SUBROUTINE tcm_define_netcdf_grid !------------------------------------------------------------------------------! ! Description: ! ------------ !> Average 3D data. !------------------------------------------------------------------------------! SUBROUTINE tcm_3d_data_averaging( mode, variable ) USE averaging, & ONLY: diss_av, kh_av, km_av USE control_parameters, & ONLY: average_count_3d IMPLICIT NONE CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< IF ( mode == 'allocate' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'diss' ) IF ( .NOT. ALLOCATED( diss_av ) ) THEN ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF diss_av = 0.0_wp CASE ( 'kh' ) IF ( .NOT. ALLOCATED( kh_av ) ) THEN ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF kh_av = 0.0_wp CASE ( 'km' ) IF ( .NOT. ALLOCATED( km_av ) ) THEN ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF km_av = 0.0_wp CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'sum' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'diss' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i) ENDDO ENDDO ENDDO CASE ( 'kh' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i) ENDDO ENDDO ENDDO CASE ( 'km' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 km_av(k,j,i) = km_av(k,j,i) + km(k,j,i) ENDDO ENDDO ENDDO CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'average' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'diss' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 diss_av(k,j,i) = diss_av(k,j,i) & / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO CASE ( 'kh' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 kh_av(k,j,i) = kh_av(k,j,i) & / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO CASE ( 'km' ) DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 km_av(k,j,i) = km_av(k,j,i) & / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO END SELECT ENDIF END SUBROUTINE tcm_3d_data_averaging !------------------------------------------------------------------------------! ! Description: ! ------------ !> Define 2D output variables. !------------------------------------------------------------------------------! SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, & two_d, nzb_do, nzt_do ) USE averaging, & ONLY: diss_av, kh_av, km_av IMPLICIT NONE CHARACTER (LEN=*) :: grid !< CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: nzb_do !< INTEGER(iwp) :: nzt_do !< LOGICAL :: found !< LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !< local !< array to which output data is resorted to found = .TRUE. SELECT CASE ( TRIM( variable ) ) CASE ( 'diss_xy', 'diss_xz', 'diss_yz' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = diss(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = diss_av(k,j,i) ENDDO ENDDO ENDDO ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'kh_xy', 'kh_xz', 'kh_yz' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = kh(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = kh_av(k,j,i) ENDDO ENDDO ENDDO ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE ( 'km_xy', 'km_xz', 'km_yz' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = km(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = km_av(k,j,i) ENDDO ENDDO ENDDO ENDIF IF ( mode == 'xy' ) grid = 'zu' CASE DEFAULT found = .FALSE. grid = 'none' END SELECT END SUBROUTINE tcm_data_output_2d !------------------------------------------------------------------------------! ! Description: ! ------------ !> Define 3D output variables. !------------------------------------------------------------------------------! SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf ) USE averaging, & ONLY: diss_av, kh_av, km_av IMPLICIT NONE CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< LOGICAL :: found !< REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !< local !< array to which output data is resorted to found = .TRUE. SELECT CASE ( TRIM( variable ) ) CASE ( 'diss' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = diss(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = diss_av(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'kh' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = kh(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = kh_av(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'km' ) IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = km(k,j,i) ENDDO ENDDO ENDDO ELSE DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = km_av(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'dummy1' ) !### remove later IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = dummy1(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'dummy2' ) !### remove later IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = dummy2(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'dummy3' ) !### remove later IF ( av == 0 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 local_pf(i,j,k) = dummy3(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE DEFAULT found = .FALSE. END SELECT END SUBROUTINE tcm_data_output_3d !------------------------------------------------------------------------------! ! Description: ! ------------ !> Allocate arrays and assign pointers. !------------------------------------------------------------------------------! SUBROUTINE tcm_init_arrays USE microphysics_mod, & ONLY: collision_turbulence USE particle_attributes, & ONLY: use_sgs_for_particles, wang_kernel IMPLICIT NONE ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) !### remove later ALLOCATE( dummy2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( dummy3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) IF ( rans_mode ) ALLOCATE( l_black(nzb:nzt+1) ) #if defined( __nopointer ) ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #else ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) #endif IF ( rans_tke_e .OR. use_sgs_for_particles .OR. wang_kernel .OR. & collision_turbulence ) THEN #if defined( __nopointer ) ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) IF ( rans_tke_e ) THEN ALLOCATE( diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF #else ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) IF ( rans_tke_e ) THEN ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ENDIF #endif ENDIF #if ! defined( __nopointer ) ! !-- Initial assignment of pointers e => e_1; e_p => e_2; te_m => e_3 IF ( rans_tke_e .OR. use_sgs_for_particles .OR. & wang_kernel .OR. collision_turbulence ) THEN diss => diss_1 IF ( rans_tke_e ) THEN diss_p => diss_2; tdiss_m => diss_3 ENDIF ENDIF #endif END SUBROUTINE tcm_init_arrays !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialization of turbulence closure module. !------------------------------------------------------------------------------! SUBROUTINE tcm_init USE arrays_3d, & ONLY: ug, vg, zu USE control_parameters, & ONLY: complex_terrain, f, kappa, dissipation_1d, topography USE model_1d_mod, & ONLY: diss1d, e1d, kh1d, km1d, l1d USE surface_mod, & ONLY: get_topography_top_index_ji IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index INTEGER(iwp) :: nz_s_shift !< INTEGER(iwp) :: nz_s_shift_l !< ! !-- Actions for initial runs IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN ! !-- Transfer initial profiles to the arrays of the 3D model DO i = nxlg, nxrg DO j = nysg, nyng e(:,j,i) = e1d kh(:,j,i) = kh1d km(:,j,i) = km1d ENDDO ENDDO IF ( constant_diffusion ) THEN e = 0.0_wp ENDIF IF ( rans_tke_e ) THEN IF ( dissipation_1d == 'prognostic' ) THEN !### Why must this be checked? DO i = nxlg, nxrg !### Should 'diss' not always DO j = nysg, nyng !### be prognostic in case rans_tke_e? diss(:,j,i) = diss1d ENDDO ENDDO ELSE DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt diss(k,j,i) = e(k,j,i) * SQRT( e(k,j,i) ) / l1d(k) ENDDO ENDDO ENDDO ENDIF ENDIF ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. & INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN IF ( constant_diffusion ) THEN km = km_constant kh = km / prandtl_number e = 0.0_wp ELSEIF ( e_init > 0.0_wp ) THEN DO k = nzb+1, nzt km(k,:,:) = c_m * l_grid(k) * SQRT( e_init ) ENDDO km(nzb,:,:) = km(nzb+1,:,:) km(nzt+1,:,:) = km(nzt,:,:) kh = km / prandtl_number e = e_init ELSE IF ( .NOT. ocean ) THEN kh = 0.01_wp ! there must exist an initial diffusion, because km = 0.01_wp ! otherwise no TKE would be produced by the ! production terms, as long as not yet ! e = (u*/cm)**2 at k=nzb+1 ELSE kh = 0.00001_wp km = 0.00001_wp ENDIF e = 0.0_wp ENDIF ENDIF ! !-- Store initial profiles for output purposes etc. hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 ) hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 ) ! !-- Initialize old and new time levels. te_m = 0.0_wp e_p = e IF ( rans_tke_e ) THEN tdiss_m = 0.0_wp diss_p = diss ENDIF ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & TRIM( initializing_actions ) == 'cyclic_fill' ) & THEN ! !-- In case of complex terrain and cyclic fill method as initialization, !-- shift initial data in the vertical direction for each point in the !-- x-y-plane depending on local surface height IF ( complex_terrain .AND. & TRIM( initializing_actions ) == 'cyclic_fill' ) THEN DO i = nxlg, nxrg DO j = nysg, nyng nz_s_shift = get_topography_top_index_ji( j, i, 's' ) e(nz_s_shift:nzt+1,j,i) = e(0:nzt+1-nz_s_shift,j,i) km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1-nz_s_shift,j,i) kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1-nz_s_shift,j,i) ENDDO ENDDO ENDIF ! !-- Initialization of the turbulence recycling method IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & turbulent_inflow ) THEN mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e ! !-- In case of complex terrain, determine vertical displacement at inflow !-- boundary and adjust mean inflow profiles IF ( complex_terrain ) THEN IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' ) ELSE nz_s_shift_l = 0 ENDIF #if defined( __parallel ) CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, & MPI_MAX, comm2d, ierr) #else nz_s_shift = nz_s_shift_l #endif mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0) ! e ENDIF ! !-- Use these mean profiles at the inflow (provided that Dirichlet !-- conditions are used) IF ( inflow_l ) THEN DO j = nysg, nyng DO k = nzb, nzt+1 e(k,j,nxlg:-1) = mean_inflow_profiles(k,5) ENDDO ENDDO ENDIF ENDIF ! !-- Inside buildings set TKE back to zero IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & topography /= 'flat' ) THEN ! !-- Inside buildings set TKE back to zero. !-- Other scalars (km, kh, diss, ...) are ignored at present, !-- maybe revise later. DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) te_m(k,j,i) = MERGE( te_m(k,j,i), 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ENDIF ! !-- Initialize new time levels (only done in order to set boundary values !-- including ghost points) e_p = e ! !-- Allthough tendency arrays are set in prognostic_equations, they have !-- have to be predefined here because they are used (but multiplied with 0) !-- there before they are set. te_m = 0.0_wp ENDIF ! !-- Calculate mixing length according to Blackadar (1962) IF ( rans_mode ) THEN IF ( f /= 0.0_wp ) THEN l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & ABS( f ) + 1E-10_wp ELSE l_max = 30.0_wp ENDIF DO k = nzb, nzt l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max ) ENDDO l_black(nzt+1) = l_black(nzt) ENDIF END SUBROUTINE tcm_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Initialize virtual velocities used later in production_e. !------------------------------------------------------------------------------! SUBROUTINE production_e_init USE arrays_3d, & ONLY: drho_air_zw, zu USE control_parameters, & ONLY: constant_flux_layer USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h IMPLICIT NONE INTEGER(iwp) :: i !< grid index x-direction INTEGER(iwp) :: j !< grid index y-direction INTEGER(iwp) :: k !< grid index z-direction INTEGER(iwp) :: m !< running index surface elements IF ( constant_flux_layer ) THEN ! !-- Calculate a virtual velocity at the surface in a way that the !-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the !-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear !-- production term at k=1 (see production_e_ij). !-- The velocity gradient has to be limited in case of too small km !-- (otherwise the timestep may be significantly reduced by large !-- surface winds). !-- not available in case of non-cyclic boundary conditions. !-- WARNING: the exact analytical solution would require the determination !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. !-- Default surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_def_h(0)%ns i = surf_def_h(0)%i(m) j = surf_def_h(0)%j(m) k = surf_def_h(0)%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_def_h(0)%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_def_h(0)%v_0(m) = v(k-1,j,i) ENDDO ! !-- Default surfaces, downward-facing surfaces !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_def_h(1)%ns i = surf_def_h(1)%i(m) j = surf_def_h(1)%j(m) k = surf_def_h(1)%k(m) surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_def_h(1)%u_0(m) = u(k+1,j,i) IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_def_h(1)%v_0(m) = v(k+1,j,i) ENDDO ! !-- Natural surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_lsm_h%ns i = surf_lsm_h%i(m) j = surf_lsm_h%j(m) k = surf_lsm_h%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_lsm_h%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_lsm_h%v_0(m) = v(k-1,j,i) ENDDO ! !-- Urban surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_usm_h%ns i = surf_usm_h%i(m) j = surf_usm_h%j(m) k = surf_usm_h%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & drho_air_zw(k-1) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_usm_h%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_usm_h%v_0(m) = v(k-1,j,i) ENDDO ENDIF END SUBROUTINE production_e_init !------------------------------------------------------------------------------! ! Description: ! ------------ !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. !> Vector-optimized version !------------------------------------------------------------------------------! SUBROUTINE tcm_prognostic USE arrays_3d, & ONLY: ddzu USE control_parameters, & ONLY: f, scalar_advec, tsc USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index INTEGER(iwp) :: m !< loop index INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp) :: sbt !< wheighting factor for sub-time step REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: advec !< advection term of TKE tendency REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: produc !< production term of TKE tendency ! !-- If required, compute prognostic equation for turbulent kinetic !-- energy (TKE) IF ( .NOT. constant_diffusion ) THEN CALL cpu_log( log_point(16), 'tke-equation', 'start' ) sbt = tsc(2) IF ( .NOT. use_upstream_for_tke ) THEN IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( e, 'e' ) ENDIF ENDIF ! !-- TKE-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN IF ( use_upstream_for_tke ) THEN tend = 0.0_wp CALL advec_s_up( e ) ELSE tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( e, 'e' ) ELSE CALL advec_s_pw( e ) ENDIF ELSE CALL advec_s_up( e ) ENDIF ENDIF ENDIF IF ( rans_tke_e ) advec = tend CALL production_e ! !-- Save production term for prognostic equation of TKE dissipation rate IF ( rans_tke_e ) produc = tend - advec IF ( .NOT. humidity ) THEN IF ( ocean ) THEN CALL diffusion_e( prho, prho_reference ) ELSE CALL diffusion_e( pt, pt_reference ) ENDIF ELSE CALL diffusion_e( vpt, pt_reference ) ENDIF ! !-- Additional sink term for flows through plant canopies IF ( plant_canopy ) CALL pcm_tendency( 6 ) CALL user_actions( 'e-tendency' ) ! !-- Prognostic equation for TKE. !-- Eliminate negative TKE values, which can occur due to numerical !-- reasons in the course of the integration. In such cases the old TKE !-- value is reduced by 90%. DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * te_m(k,j,i) ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) ENDDO ENDDO ENDDO ! !-- Use special boundary condition in case of TKE-e closure IF ( rans_tke_e ) THEN DO i = nxl, nxr DO j = nys, nyn surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2 ENDDO ENDDO ENDDO ENDIF ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt te_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt te_m(k,j,i) = -9.5625_wp * tend(k,j,i) & + 5.3125_wp * te_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(16), 'tke-equation', 'stop' ) ENDIF ! TKE equation ! !-- If required, compute prognostic equation for TKE dissipation rate IF ( rans_tke_e ) THEN CALL cpu_log( log_point(33), 'diss-equation', 'start' ) sbt = tsc(2) IF ( .NOT. use_upstream_for_tke ) THEN IF ( scalar_advec == 'bc-scheme' ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN ! !-- Bott-Chlond scheme always uses Euler time step. Thus: sbt = 1.0_wp ENDIF tend = 0.0_wp CALL advec_s_bc( diss, 'diss' ) ENDIF ENDIF ! !-- dissipation-tendency terms with no communication IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN IF ( use_upstream_for_tke ) THEN tend = 0.0_wp CALL advec_s_up( diss ) ELSE tend = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( diss, 'diss' ) ELSE CALL advec_s_pw( diss ) ENDIF ELSE CALL advec_s_up( diss ) ENDIF ENDIF ENDIF ! !-- Production of TKE dissipation rate DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k) tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h & !### needs revision / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) & * SQRT(e(k,j,i)) * produc(k,j,i) ENDDO ENDDO ENDDO CALL diffusion_diss ! !-- Additional sink term for flows through plant canopies ! IF ( plant_canopy ) CALL pcm_tendency( ? ) !### what to do with this? ! CALL user_actions( 'diss-tendency' ) !### not yet implemented ! !-- Prognostic equation for TKE dissipation. !-- Eliminate negative dissipation values, which can occur due to numerical !-- reasons in the course of the integration. In such cases the old !-- dissipation value is reduced by 90%. DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & tsc(3) * tdiss_m(k,j,i) ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) IF ( diss_p(k,j,i) < 0.0_wp ) & diss_p(k,j,i) = 0.1_wp * diss(k,j,i) ENDDO ENDDO ENDDO ! !-- Use special boundary condition in case of TKE-e closure DO i = nxl, nxr DO j = nys, nyn surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / kappa * ddzu(k) ENDDO ENDDO ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tdiss_m(k,j,i) = tend(k,j,i) ENDDO ENDDO ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) & + 5.3125_wp * tdiss_m(k,j,i) ENDDO ENDDO ENDDO ENDIF ENDIF CALL cpu_log( log_point(33), 'diss-equation', 'stop' ) ENDIF END SUBROUTINE tcm_prognostic !------------------------------------------------------------------------------! ! Description: ! ------------ !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. !> Cache-optimized version !------------------------------------------------------------------------------! SUBROUTINE tcm_prognostic_ij( i, j, i_omp, tn ) USE arrays_3d, & ONLY: ddzu, diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, & flux_l_diss, flux_l_e, flux_s_diss, flux_s_e USE control_parameters, & ONLY: f, tsc USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< loop index x direction INTEGER(iwp) :: i_omp !< INTEGER(iwp) :: j !< loop index y direction INTEGER(iwp) :: k !< loop index z direction INTEGER(iwp) :: m !< loop index INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position INTEGER(iwp) :: tn !< REAL(wp), DIMENSION(nzb:nzt+1) :: advec !< advection term of TKE tendency REAL(wp), DIMENSION(nzb:nzt+1) :: produc !< production term of TKE tendency ! !-- If required, compute prognostic equation for turbulent kinetic !-- energy (TKE) IF ( .NOT. constant_diffusion ) THEN ! !-- Tendency-terms for TKE tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' & .AND. .NOT. use_upstream_for_tke ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, & flux_l_e, diss_l_e , i_omp, tn ) ELSE CALL advec_s_pw( i, j, e ) ENDIF ELSE CALL advec_s_up( i, j, e ) ENDIF advec(:) = tend(:,j,i) CALL production_e( i, j ) produc(:) = tend(:,j,i) - advec(:) IF ( .NOT. humidity ) THEN IF ( ocean ) THEN CALL diffusion_e( i, j, prho, prho_reference ) ELSE CALL diffusion_e( i, j, pt, pt_reference ) ENDIF ELSE CALL diffusion_e( i, j, vpt, pt_reference ) ENDIF ! !-- Additional sink term for flows through plant canopies IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 ) CALL user_actions( i, j, 'e-tendency' ) ! !-- Prognostic equation for TKE. !-- Eliminate negative TKE values, which can occur due to numerical !-- reasons in the course of the integration. In such cases the old !-- TKE value is reduced by 90%. DO k = nzb+1, nzt e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & tsc(3) * te_m(k,j,i) ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) IF ( e_p(k,j,i) <= 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) ENDDO ! !-- Use special boundary condition in case of TKE-e closure IF ( rans_tke_e ) THEN surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2 ENDDO ENDIF ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt te_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 5.3125_wp * te_m(k,j,i) ENDDO ENDIF ENDIF ENDIF ! TKE equation ! !-- If required, compute prognostic equation for TKE dissipation rate IF ( rans_tke_e ) THEN ! !-- Tendency-terms for dissipation tend(:,j,i) = 0.0_wp IF ( timestep_scheme(1:5) == 'runge' & .AND. .NOT. use_upstream_for_tke ) THEN IF ( ws_scheme_sca ) THEN CALL advec_s_ws( i, j, diss, 'diss', flux_s_diss, diss_s_diss, & flux_l_diss, diss_l_diss, i_omp, tn ) ELSE CALL advec_s_pw( i, j, diss ) ENDIF ELSE CALL advec_s_up( i, j, diss ) ENDIF ! !-- Production of TKE dissipation rate DO k = nzb+1, nzt ! tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k) tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h & !### needs revision / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) & * SQRT(e(k,j,i)) * produc(k) ENDDO CALL diffusion_diss( i, j ) ! !-- Additional sink term for flows through plant canopies ! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !### not yet implemented ! CALL user_actions( i, j, 'diss-tendency' ) !### not yet implemented ! !-- Prognostic equation for TKE dissipation !-- Eliminate negative dissipation values, which can occur due to !-- numerical reasons in the course of the integration. In such cases !-- the old dissipation value is reduced by 90%. DO k = nzb+1, nzt diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & tsc(3) * tdiss_m(k,j,i) ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 )& ) IF ( diss_p(k,j,i) <= 0.0_wp ) diss_p(k,j,i) = 0.1_wp * diss(k,j,i) ENDDO ! !-- Use special boundary condition in case of TKE-e closure surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / kappa * ddzu(k) ENDDO ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt tdiss_m(k,j,i) = tend(k,j,i) ENDDO ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 5.3125_wp * tdiss_m(k,j,i) ENDDO ENDIF ENDIF ! IF ( intermediate_timestep_count == 1 ) dummy1(:,j,i) = e_p(:,j,i) ! IF ( intermediate_timestep_count == 2 ) dummy2(:,j,i) = e_p(:,j,i) ! IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = e_p(:,j,i) ENDIF ! dissipation equation END SUBROUTINE tcm_prognostic_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Production terms (shear + buoyancy) of the TKE. !> Vector-optimized version !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is !> not considered well! !------------------------------------------------------------------------------! SUBROUTINE production_e USE arrays_3d, & ONLY: ddzw, dd2zu, drho_air_zw, q, ql USE cloud_parameters, & ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt USE control_parameters, & ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, neutral, & rho_reference, use_single_reference_value, use_surface_fluxes, & use_top_fluxes USE grid_variables, & ONLY: ddx, dx, ddy, dy USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< running index x-direction INTEGER(iwp) :: j !< running index y-direction INTEGER(iwp) :: k !< running index z-direction INTEGER(iwp) :: l !< running index for different surface type orientation INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp) :: def !< REAL(wp) :: flag !< flag to mask topography REAL(wp) :: k1 !< REAL(wp) :: k2 !< REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces REAL(wp) :: theta !< REAL(wp) :: temp !< REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation REAL(wp) :: usvs !< momentum flux u"v" REAL(wp) :: vsus !< momentum flux v"u" REAL(wp) :: wsus !< momentum flux w"u" REAL(wp) :: wsvs !< momentum flux w"v" REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudx !< Gradient of u-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudy !< Gradient of u-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudz !< Gradient of u-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdx !< Gradient of v-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdy !< Gradient of v-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdz !< Gradient of v-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdx !< Gradient of w-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdy !< Gradient of w-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdz !< Gradient of w-component in z-direction DO i = nxl, nxr IF ( constant_flux_layer ) THEN ! !-- Calculate TKE production by shear. Calculate gradients at all grid !-- points first, gradients at surface-bounded grid points will be !-- overwritten further below. DO j = nys, nyn DO k = nzb+1, nzt dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * & dd2zu(k) dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * & dd2zu(k) dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ENDDO ENDDO ! !-- Position beneath wall !-- (2) - Will allways be executed. !-- 'bottom and wall: use u_0,v_0 and wall functions' DO j = nys, nyn ! !-- Compute gradients at north- and south-facing surfaces. !-- First, for default surfaces, then for urban surfaces. !-- Note, so far no natural vertical surfaces implemented DO l = 0, 1 surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) usvs = surf_def_v(l)%mom_flux_tke(0,m) wsvs = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) usvs = surf_lsm_v(l)%mom_flux_tke(0,m) wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) usvs = surf_usm_v(l)%mom_flux_tke(0,m) wsvs = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at east- and west-facing walls DO l = 2, 3 surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) vsus = surf_def_v(l)%mom_flux_tke(0,m) wsus = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) vsus = surf_lsm_v(l)%mom_flux_tke(0,m) wsus = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) vsus = surf_usm_v(l)%mom_flux_tke(0,m) wsus = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at upward-facing surfaces surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Compute gradients at downward-facing walls, only for !-- non-natural default surfaces surf_s = surf_def_h(1)%start_index(j,i) surf_e = surf_def_h(1)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(1)%k(m) dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) ENDDO ENDDO DO j = nys, nyn DO k = nzb+1, nzt def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + & dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & dwdy(k,j)*dvdz(k,j) ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDDO ELSE DO j = nys, nyn ! !-- Calculate TKE production by shear. Here, no additional !-- wall-bounded code is considered. !-- Why? DO k = nzb+1, nzt dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * & dd2zu(k) dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * & dd2zu(k) dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * & ddzw(k) def = 2.0_wp * ( & dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 & ) + & dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2.0_wp * ( & dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & dwdy(k,j)*dvdz(k,j) & ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 29 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDDO ENDIF ! !-- If required, calculate TKE production by buoyancy IF ( .NOT. neutral ) THEN IF ( .NOT. humidity ) THEN IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux !-- in the bottom and top surface layer DO j = nys, nyn DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO ! !-- Treatment of near-surface grid points, at up- and down- !-- ward facing surfaces IF ( use_surface_fluxes ) THEN DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & drho_air_zw(k-1) * & surf_def_h(l)%shf(m) ENDDO ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & drho_air_zw(k) * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDDO ELSE DO j = nys, nyn DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) & * drho_air_zw(k-1) & * surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) & * drho_air_zw(k-1) & * surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) & * drho_air_zw(k-1) & * surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) & * drho_air_zw(k) & * surf_def_h(2)%shf(m) ENDDO ENDIF ENDDO ENDIF ELSE DO j = nys, nyn DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & pt(k,j,i) * ( ql(k+1,j,i) - & ql(k-1,j,i) ) ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDIF ENDDO ENDDO IF ( use_surface_fluxes ) THEN DO j = nys, nyn ! !-- Treat horizontal default surfaces DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_def_h(l)%shf(m) + & k2 * surf_def_h(l)%qsws(m) & ) * drho_air_zw(k-1) ENDDO ENDDO ! !-- Treat horizontal natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_lsm_h%shf(m) + & k2 * surf_lsm_h%qsws(m) & ) * drho_air_zw(k-1) ENDDO ! !-- Treat horizontal urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_usm_h%shf(m) + & k2 * surf_usm_h%qsws(m) & ) * drho_air_zw(k-1) ENDDO ENDDO ENDIF IF ( use_top_fluxes ) THEN DO j = nys, nyn surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_def_h(2)%shf(m) + & k2 * surf_def_h(2)%qsws(m) & ) * drho_air_zw(k) ENDDO ENDDO ENDIF ENDIF ENDIF ENDDO END SUBROUTINE production_e !------------------------------------------------------------------------------! ! Description: ! ------------ !> Production terms (shear + buoyancy) of the TKE. !> Cache-optimized version !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is !> not considered well! !------------------------------------------------------------------------------! SUBROUTINE production_e_ij( i, j ) USE arrays_3d, & ONLY: ddzw, dd2zu, drho_air_zw, q, ql USE cloud_parameters, & ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt USE control_parameters, & ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, neutral, & rho_reference, use_single_reference_value, use_surface_fluxes, & use_top_fluxes USE grid_variables, & ONLY: ddx, dx, ddy, dy USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< running index x-direction INTEGER(iwp) :: j !< running index y-direction INTEGER(iwp) :: k !< running index z-direction INTEGER(iwp) :: l !< running index for different surface type orientation INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp) :: def !< REAL(wp) :: flag !< flag to mask topography REAL(wp) :: k1 !< REAL(wp) :: k2 !< REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces REAL(wp) :: theta !< REAL(wp) :: temp !< REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation REAL(wp) :: usvs !< momentum flux u"v" REAL(wp) :: vsus !< momentum flux v"u" REAL(wp) :: wsus !< momentum flux w"u" REAL(wp) :: wsvs !< momentum flux w"v" REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction IF ( constant_flux_layer ) THEN ! !-- Calculate TKE production by shear. Calculate gradients at all grid !-- points first, gradients at surface-bounded grid points will be !-- overwritten further below. DO k = nzb+1, nzt dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ENDDO ! !-- Compute gradients at north- and south-facing surfaces. !-- Note, no vertical natural surfaces so far. DO l = 0, 1 ! !-- Default surfaces surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) usvs = surf_def_v(l)%mom_flux_tke(0,m) wsvs = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) usvs = surf_lsm_v(l)%mom_flux_tke(0,m) wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) usvs = surf_usm_v(l)%mom_flux_tke(0,m) wsvs = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at east- and west-facing walls DO l = 2, 3 ! !-- Default surfaces surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) vsus = surf_def_v(l)%mom_flux_tke(0,m) wsus = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) vsus = surf_lsm_v(l)%mom_flux_tke(0,m) wsus = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) vsus = surf_usm_v(l)%mom_flux_tke(0,m) wsus = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at upward-facing walls, first for !-- non-natural default surfaces surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Compute gradients at downward-facing walls, only for !-- non-natural default surfaces surf_s = surf_def_h(1)%start_index(j,i) surf_e = surf_def_h(1)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(1)%k(m) dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) ENDDO DO k = nzb+1, nzt def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) ! !-- Production term according to Kato and Launder (1993) ! def = SQRT( ( dudx(k) + dudy(k) + dudz(k) + & ! dvdx(k) + dvdy(k) + dvdz(k) + & ! dwdx(k) + dwdy(k) + dwdz(k) & ! )**4 - & ! ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 + & ! 2.0_wp * ( dudy(k) * dvdx(k) + & ! dudz(k) * dwdx(k) + & ! dvdz(k) * dwdy(k) ) & ! )**2 & ! ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ELSE ! !-- Calculate TKE production by shear. Here, no additional !-- wall-bounded code is considered. !-- Why? DO k = nzb+1, nzt dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) ! !-- Production term according to Kato and Launder (1993) ! def = SQRT( ( dudx(k) + dudy(k) + dudz(k) + & ! dvdx(k) + dvdy(k) + dvdz(k) + & ! dwdx(k) + dwdy(k) + dwdz(k) & ! )**4 - & ! ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 + & ! 2.0_wp * ( dudy(k) * dvdx(k) + & ! dudz(k) * dwdx(k) + & ! dvdz(k) * dwdy(k) ) & ! )**2 & ! ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 29 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDIF ! !-- If required, calculate TKE production by buoyancy IF ( .NOT. neutral ) THEN IF ( .NOT. humidity ) THEN IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux in !-- the bottom and top surface layer DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & drho_air_zw(k-1) * & surf_def_h(l)%shf(m) ENDDO ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( rho_reference, prho(k,j,i), & use_single_reference_value ) * & drho_air_zw(k) * & surf_def_h(2)%shf(m) ENDDO ENDIF ELSE DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & drho_air_zw(k-1) * & surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & drho_air_zw(k-1) * & surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & drho_air_zw(k-1) * & surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / & MERGE( pt_reference, pt(k,j,i), & use_single_reference_value ) * & drho_air_zw(k) * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDIF ELSE DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & pt(k,j,i) * ( ql(k+1,j,i) - & ql(k-1,j,i) ) ) * dd2zu(k) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDIF ENDDO IF ( use_surface_fluxes ) THEN ! !-- Treat horizontal default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_def_h(l)%shf(m) + & k2 * surf_def_h(l)%qsws(m) & ) * drho_air_zw(k-1) ENDDO ENDDO ! !-- Treat horizontal natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_lsm_h%shf(m) + & k2 * surf_lsm_h%qsws(m) & ) * drho_air_zw(k-1) ENDDO ! !-- Treat horizontal urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1 * surf_usm_h%shf(m) + & k2 * surf_usm_h%qsws(m) & ) * drho_air_zw(k-1) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / & MERGE( vpt_reference, vpt(k,j,i), & use_single_reference_value ) * & ( k1* surf_def_h(2)%shf(m) + & k2 * surf_def_h(2)%qsws(m) & ) * drho_air_zw(k) ENDDO ENDIF ENDIF ENDIF END SUBROUTINE production_e_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Diffusion and dissipation terms for the TKE. !> Vector-optimized version !------------------------------------------------------------------------------! SUBROUTINE diffusion_e( var, var_reference ) USE arrays_3d, & ONLY: ddzu, ddzw, drho_air, rho_air_zw USE grid_variables, & ONLY: ddx2, ddy2 USE microphysics_mod, & ONLY: collision_turbulence USE particle_attributes, & ONLY: use_sgs_for_particles, wang_kernel USE surface_mod, & ONLY : bc_h IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction INTEGER(iwp) :: m !< running index surface elements REAL(wp) :: flag !< flag to mask topography REAL(wp) :: l !< mixing length REAL(wp) :: ll !< adjusted l REAL(wp) :: var_reference !< #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< #else REAL(wp), DIMENSION(:,:,:), POINTER :: var !< #endif REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !< TKE dissipation ! !-- Calculate the tendency terms DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Calculate dissipation IF ( les_mw ) THEN CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) dissipation(k,j) = ( 0.19_wp + 0.74_wp * l / ll ) & * e(k,j,i) * SQRT( e(k,j,i) ) / l ELSEIF ( rans_tke_l ) THEN CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) dissipation(k,j) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll ELSEIF ( rans_tke_e ) THEN dissipation(k,j) = diss(k,j,i) ENDIF ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) & + ( & ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & ) * ddx2 * flag & + ( & ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & ) * ddy2 * flag & + ( & ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & * rho_air_zw(k) & - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & * rho_air_zw(k-1) & ) * ddzw(k) * drho_air(k) * flag & - dissipation(k,j) * flag ENDDO ENDDO ! !-- Store dissipation if needed for calculating the sgs particle !-- velocities IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. & wang_kernel .OR. collision_turbulence ) ) THEN DO j = nys, nyn DO k = nzb+1, nzt diss(k,j,i) = dissipation(k,j) * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDIF ENDDO ! !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. & wang_kernel .OR. collision_turbulence ) ) THEN ! !-- Upward facing surfaces DO m = 1, bc_h(0)%ns i = bc_h(0)%i(m) j = bc_h(0)%j(m) k = bc_h(0)%k(m) diss(k-1,j,i) = diss(k,j,i) ENDDO ! !-- Downward facing surfaces DO m = 1, bc_h(1)%ns i = bc_h(1)%i(m) j = bc_h(1)%j(m) k = bc_h(1)%k(m) diss(k+1,j,i) = diss(k,j,i) ENDDO ENDIF END SUBROUTINE diffusion_e !------------------------------------------------------------------------------! ! Description: ! ------------ !> Diffusion and dissipation terms for the TKE. !> Cache-optimized version !------------------------------------------------------------------------------! SUBROUTINE diffusion_e_ij( i, j, var, var_reference ) USE arrays_3d, & ONLY: ddzu, ddzw, drho_air, rho_air_zw USE grid_variables, & ONLY: ddx2, ddy2 USE microphysics_mod, & ONLY: collision_turbulence USE particle_attributes, & ONLY: use_sgs_for_particles, wang_kernel USE surface_mod, & ONLY : bc_h IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint REAL(wp) :: flag !< flag to mask topography REAL(wp) :: l !< mixing length REAL(wp) :: ll !< adjusted l REAL(wp) :: var_reference !< #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< #else REAL(wp), DIMENSION(:,:,:), POINTER :: var !< #endif REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !< dissipation of TKE ! !-- Calculate the mixing length (for dissipation) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Calculate dissipation... !-- ...in case of LES IF ( les_mw ) THEN CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) dissipation(k) = ( 0.19_wp + 0.74_wp * l / ll ) & * e(k,j,i) * SQRT( e(k,j,i) ) / l ! !-- ...in case of RANS ELSEIF ( rans_tke_l ) THEN CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) dissipation(k) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll ELSEIF ( rans_tke_e ) THEN dissipation(k) = diss(k,j,i) ENDIF ! !-- Calculate the tendency term tend(k,j,i) = tend(k,j,i) & + ( & ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & ) * ddx2 * flag / sig_e & + ( & ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & ) * ddy2 * flag / sig_e & + ( & ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & * rho_air_zw(k) & - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & * rho_air_zw(k-1) & ) * ddzw(k) * drho_air(k) * flag / sig_e & - dissipation(k) * flag ENDDO ! !-- Store dissipation if needed for calculating the sgs particle velocities IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. wang_kernel & .OR. collision_turbulence ) ) THEN DO k = nzb+1, nzt diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ! !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) !-- For each surface type determine start and end index (in case of elevated !-- topography several up/downward facing surfaces may exist. surf_s = bc_h(0)%start_index(j,i) surf_e = bc_h(0)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(0)%k(m) diss(k-1,j,i) = diss(k,j,i) ENDDO ! !-- Downward facing surfaces surf_s = bc_h(1)%start_index(j,i) surf_e = bc_h(1)%end_index(j,i) DO m = surf_s, surf_e k = bc_h(1)%k(m) diss(k+1,j,i) = diss(k,j,i) ENDDO ENDIF END SUBROUTINE diffusion_e_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Diffusion term for the TKE dissipation rate !> Vector-optimized version !------------------------------------------------------------------------------! SUBROUTINE diffusion_diss() USE arrays_3d, & ONLY: ddzu, ddzw, drho_air, rho_air_zw USE grid_variables, & ONLY: ddx2, ddy2 IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction REAL(wp) :: flag !< flag to mask topography ! !-- Calculate the tendency terms DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) & + ( & ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) & - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) & ) * ddx2 * flag & + ( & ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) & - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) & ) * ddy2 * flag & + ( & ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) & * rho_air_zw(k) & - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) & * rho_air_zw(k-1) & ) * ddzw(k) * drho_air(k) * flag & - c_2 * diss(k,j,i)**2 & / ( e(k,j,i) + 1.0E-20_wp ) * flag ENDDO ENDDO ENDDO END SUBROUTINE diffusion_diss !------------------------------------------------------------------------------! ! Description: ! ------------ !> Diffusion term for the TKE dissipation rate !> Cache-optimized version !------------------------------------------------------------------------------! SUBROUTINE diffusion_diss_ij( i, j ) USE arrays_3d, & ONLY: ddzu, ddzw, drho_air, rho_air_zw USE grid_variables, & ONLY: ddx2, ddy2 IMPLICIT NONE INTEGER(iwp) :: i !< running index x direction INTEGER(iwp) :: j !< running index y direction INTEGER(iwp) :: k !< running index z direction REAL(wp) :: flag !< flag to mask topography REAL(wp), DIMENSION(nzb+1:nzt) :: tend_temp ! !-- Calculate the mixing length (for dissipation) DO k = nzb+1, nzt ! !-- Predetermine flag to mask topography flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Calculate the tendency term tend_temp(k) = ( & ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) & - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) & ) * ddx2 * flag / sig_diss & + ( & ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) & - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) & ) * ddy2 * flag / sig_diss & + ( & ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) & * rho_air_zw(k) & - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) & * rho_air_zw(k-1) & ) * ddzw(k) * drho_air(k) * flag / sig_diss & - c_2 * diss(k,j,i)**2 & / ( e(k,j,i) + 1.0E-20_wp ) * flag tend(k,j,i) = tend(k,j,i) + tend_temp(k) ENDDO END SUBROUTINE diffusion_diss_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate mixing length for LES mode. !------------------------------------------------------------------------------! SUBROUTINE mixing_length_les( i, j, k, l, ll, var, var_reference ) USE arrays_3d, & ONLY: dd2zu, l_grid, l_wall USE control_parameters, & ONLY: atmos_ocean_sign, g, use_single_reference_value, & wall_adjustment, wall_adjustment_factor IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: dvar_dz !< vertical gradient of var REAL(wp) :: l !< mixing length REAL(wp) :: l_stable !< mixing length according to stratification REAL(wp) :: ll !< adjusted l_grid REAL(wp) :: var_reference !< var at reference height #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature #else REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature #endif dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) IF ( dvar_dz > 0.0_wp ) THEN IF ( use_single_reference_value ) THEN l_stable = 0.76_wp * SQRT( e(k,j,i) ) & / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp ELSE l_stable = 0.76_wp * SQRT( e(k,j,i) ) & / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp ENDIF ELSE l_stable = l_grid(k) ENDIF ! !-- Adjustment of the mixing length IF ( wall_adjustment ) THEN l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable ) ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) ) ELSE l = MIN( l_grid(k), l_stable ) ll = l_grid(k) ENDIF END SUBROUTINE mixing_length_les !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate mixing length for RANS mode. !------------------------------------------------------------------------------! SUBROUTINE mixing_length_rans( i, j, k, l, l_diss, var, var_reference ) USE arrays_3d, & ONLY: dd2zu USE control_parameters, & ONLY: atmos_ocean_sign, g, use_single_reference_value IMPLICIT NONE INTEGER(iwp) :: i !< loop index INTEGER(iwp) :: j !< loop index INTEGER(iwp) :: k !< loop index REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector REAL(wp) :: dvar_dz !< vertical gradient of var REAL(wp) :: l !< mixing length REAL(wp) :: l_diss !< mixing length for dissipation REAL(wp) :: rif !< Richardson flux number REAL(wp) :: var_reference !< var at reference height #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature #else REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature #endif dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) duv2_dz2 = ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2 & + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 & + 1E-30_wp IF ( use_single_reference_value ) THEN rif = g / var_reference * dvar_dz / duv2_dz2 ELSE rif = g / var(k,j,i) * dvar_dz / duv2_dz2 ENDIF rif = MAX( rif, -5.0_wp ) rif = MIN( rif, 1.0_wp ) ! !-- Calculate diabatic mixing length using Dyer-profile functions IF ( rif >= 0.0_wp ) THEN l = l_black(k) / ( 1.0_wp + 5.0_wp * rif ) l_diss = l ELSE ! !-- In case of unstable stratification, use mixing length of neutral case !-- for l, but consider profile functions for l_diss l = l_black(k) l_diss = l_black(k) * SQRT( 1.0_wp - 16.0_wp * rif ) ENDIF END SUBROUTINE mixing_length_rans !------------------------------------------------------------------------------! ! Description: ! ------------ !> Computation of the turbulent diffusion coefficients for momentum and heat !> according to Prandtl-Kolmogorov. !------------------------------------------------------------------------------! SUBROUTINE tcm_diffusivities( var, var_reference ) USE control_parameters, & ONLY: e_min, outflow_l, outflow_n, outflow_r, outflow_s USE statistics, & ONLY : rmask, sums_l_l USE surface_mod, & ONLY : bc_h, surf_def_h IMPLICIT NONE INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< INTEGER(iwp) :: m !< INTEGER(iwp) :: n !< INTEGER(iwp) :: omp_get_thread_num !< INTEGER(iwp) :: sr !< INTEGER(iwp) :: tn !< REAL(wp) :: flag !< REAL(wp) :: l !< REAL(wp) :: ll !< REAL(wp) :: var_reference !< #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< #else REAL(wp), DIMENSION(:,:,:), POINTER :: var !< #endif ! !-- Default thread number in case of one thread tn = 0 ! !-- Initialization for calculation of the mixing length profile sums_l_l = 0.0_wp ! !-- Compute the turbulent diffusion coefficient for momentum !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag) !$ tn = omp_get_thread_num() ! !-- Introduce an optional minimum tke IF ( e_min > 0.0_wp ) THEN !$OMP DO DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt e(k,j,i) = MAX( e(k,j,i), e_min ) * & MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ENDDO ENDDO ENDDO ENDIF IF ( les_mw ) THEN !$OMP DO DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Determine the mixing length for LES closure CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) ! !-- Compute diffusion coefficients for momentum and heat km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag ! !-- Summation for averaged profile (cf. flow_statistics) DO sr = 0, statistic_regions sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & * flag ENDDO ENDDO ENDDO ENDDO ELSEIF ( rans_tke_l ) THEN !$OMP DO DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Mixing length for RANS mode with TKE-l closure CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) ! !-- Compute diffusion coefficients for momentum and heat km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag kh(k,j,i) = km(k,j,i) / prandtl_number * flag ! !-- Summation for averaged profile (cf. flow_statistics) DO sr = 0, statistic_regions sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & * flag ENDDO ENDDO ENDDO ENDDO ELSEIF ( rans_tke_e ) THEN !$OMP DO DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb+1, nzt flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) ! !-- Compute diffusion coefficients for momentum and heat km(k,j,i) = c_mu * e(k,j,i)**2 / ( diss(k,j,i) + 1.E-10 ) * flag kh(k,j,i) = km(k,j,i) / prandtl_number * flag ENDDO ENDDO ENDDO ENDIF sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for ! data output !$OMP END PARALLEL ! !-- Set vertical boundary values (Neumann conditions both at upward- and !-- downward facing walls. To set wall-boundary values, the surface data type !-- is applied. !-- Horizontal boundary conditions at vertical walls are not set because !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary !-- values of the diffusivities are not needed. IF ( .NOT. rans_tke_e ) THEN ! !-- Upward-facing !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(0)%ns i = bc_h(0)%i(m) j = bc_h(0)%j(m) k = bc_h(0)%k(m) km(k-1,j,i) = km(k,j,i) kh(k-1,j,i) = kh(k,j,i) ENDDO ! !-- Downward facing surfaces !$OMP PARALLEL DO PRIVATE( i, j, k ) DO m = 1, bc_h(1)%ns i = bc_h(1)%i(m) j = bc_h(1)%j(m) k = bc_h(1)%k(m) km(k+1,j,i) = km(k,j,i) kh(k+1,j,i) = kh(k,j,i) ENDDO ELSE ! !-- Up- and downward facing surfaces DO n = 0, 1 DO m = 1, surf_def_h(n)%ns i = surf_def_h(n)%i(m) j = surf_def_h(n)%j(m) k = surf_def_h(n)%k(m) km(k,j,i) = kappa * surf_def_h(n)%us(m) * dzu(k) kh(k,j,i) = 1.35_wp * km(k,j,i) ENDDO ENDDO CALL exchange_horiz( km, nbgp ) CALL exchange_horiz( kh, nbgp ) ENDIF ! !-- Model top !$OMP PARALLEL DO DO i = nxlg, nxrg DO j = nysg, nyng km(nzt+1,j,i) = km(nzt,j,i) kh(nzt+1,j,i) = kh(nzt,j,i) ENDDO ENDDO ! !-- Set Neumann boundary conditions at the outflow boundaries in case of !-- non-cyclic lateral boundaries IF ( outflow_l ) THEN km(:,:,nxl-1) = km(:,:,nxl) kh(:,:,nxl-1) = kh(:,:,nxl) ENDIF IF ( outflow_r ) THEN km(:,:,nxr+1) = km(:,:,nxr) kh(:,:,nxr+1) = kh(:,:,nxr) ENDIF IF ( outflow_s ) THEN km(:,nys-1,:) = km(:,nys,:) kh(:,nys-1,:) = kh(:,nys,:) ENDIF IF ( outflow_n ) THEN km(:,nyn+1,:) = km(:,nyn,:) kh(:,nyn+1,:) = kh(:,nyn,:) ENDIF END SUBROUTINE tcm_diffusivities !------------------------------------------------------------------------------! ! Description: ! ------------ !> Swapping of timelevels. !------------------------------------------------------------------------------! SUBROUTINE tcm_swap_timelevel ( mod_count ) IMPLICIT NONE INTEGER(iwp) :: i !< loop index x direction INTEGER(iwp) :: j !< loop index y direction INTEGER(iwp) :: k !< loop index z direction INTEGER, INTENT(IN) :: mod_count !< #if defined( __nopointer ) IF ( .NOT. constant_diffusion ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 e(k,j,i) = e_p(k,j,i) ENDDO ENDDO ENDDO ENDIF IF ( rans_tke_e ) THEN DO i = nxlg, nxrg DO j = nysg, nyng DO k = nzb, nzt+1 diss(k,j,i) = diss_p(k,j,i) ENDDO ENDDO ENDDO ENDIF #else SELECT CASE ( mod_count ) CASE ( 0 ) IF ( .NOT. constant_diffusion ) THEN e => e_1; e_p => e_2 ENDIF IF ( rans_tke_e ) THEN diss => diss_1; diss_p => diss_2 ENDIF CASE ( 1 ) IF ( .NOT. constant_diffusion ) THEN e => e_2; e_p => e_1 ENDIF IF ( rans_tke_e ) THEN diss => diss_2; diss_p => diss_1 ENDIF END SELECT #endif END SUBROUTINE tcm_swap_timelevel END MODULE turbulence_closure_mod