Changeset 3130


Ignore:
Timestamp:
Jul 16, 2018 11:08:55 AM (3 years ago)
Author:
gronemeier
Message:

calculate km according to MOST within the surface layer

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r3129 r3130  
    2525# -----------------
    2626# $Id$
     27# add surface_layer_fluxes_mod to turbulence_closure_mod
     28#
     29# 3129 2018-07-16 07:45:13Z gronemeier
    2730# add turbulence_closure_mod to parin
    2831#
     
    15311534        plant_canopy_model_mod.o \
    15321535        pmc_interface_mod.o \
     1536        surface_layer_fluxes_mod.o \
    15331537        surface_mod.o \
    15341538        user_actions.o
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r3045 r3130  
    2626! -----------------
    2727! $Id$
     28! move phi_m from turbulence_closure_mod
     29!
     30! 3045 2018-05-28 07:55:41Z Giersch
    2831! Error message revised
    2932!
     
    287290    PRIVATE
    288291
    289     PUBLIC init_surface_layer_fluxes, surface_layer_fluxes
     292    PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
    290293
    291294    INTERFACE init_surface_layer_fluxes
    292295       MODULE PROCEDURE init_surface_layer_fluxes
    293296    END INTERFACE init_surface_layer_fluxes
     297
     298    INTERFACE phi_m
     299       MODULE PROCEDURE phi_m
     300    END INTERFACE phi_m
    294301
    295302    INTERFACE surface_layer_fluxes
     
    22932300    END FUNCTION psi_h
    22942301
     2302
     2303!------------------------------------------------------------------------------!
     2304! Description:
     2305! ------------
     2306!> Calculates stability function for momentum
     2307!>
     2308!> @author Hauke Wurps
     2309!------------------------------------------------------------------------------!
     2310    FUNCTION phi_m( zeta )
     2311   
     2312       IMPLICIT NONE
     2313   
     2314       REAL(wp)            :: phi_m         !< Value of the function
     2315       REAL(wp)            :: zeta          !< Stability parameter z/L
     2316   
     2317       REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
     2318       REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
     2319   
     2320       IF ( zeta < 0.0_wp )  THEN
     2321          phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
     2322       ELSE
     2323          phi_m = 1.0_wp + c * zeta
     2324       ENDIF
     2325   
     2326    END FUNCTION phi_m
     2327
    22952328 END MODULE surface_layer_fluxes_mod
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3129 r3130  
    2525! -----------------
    2626! $Id$
     27! - move boundary condition of km and kh to tcm_diffusivities
     28! - calculate km at boundaries according to MOST
     29! - move phi_m to surface_layer_fluxes_mod
     30!
     31! 3129 2018-07-16 07:45:13Z gronemeier
    2732! - move limitation of diss to boundary_conds
    2833! - move boundary conditions for e and diss to boundary_conds
     
    19101915    INTEGER(iwp) ::  m      !< running index surface elements
    19111916
    1912     REAL(wp)     :: km_sfc  !< km according to MOST
    1913 
    19141917    IF ( constant_flux_layer )  THEN
    19151918!
     
    19221925!--    surface winds).
    19231926!--    not available in case of non-cyclic boundary conditions.
    1924 !--    WARNING: the exact analytical solution would require the determination
    1925 !--             of the eddy diffusivity by km = u* * kappa * zp / phi_m.
    19261927!--    Default surfaces, upward-facing
    19271928       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     
    19321933          k = surf_def_h(0)%k(m)
    19331934!
    1934 !--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
     1935!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    19351936!--       and km are not on the same grid. Actually, a further
    19361937!--       interpolation of km onto the u/v-grid is necessary. However, the
    19371938!--       effect of this error is negligible.
    1938           km_sfc = surf_def_h(0)%us(m) * kappa * surf_def_h(0)%z_mo(m) /       &
    1939                    phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
    1940 
    19411939          surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
    19421940                                     drho_air_zw(k-1)               *          &
    19431941                                     ( zu(k+1) - zu(k-1)    )       /          &
    1944                                      ( km_sfc  + 1.0E-20_wp )
     1942                                     ( km(k,j,i)  + 1.0E-20_wp )
    19451943          surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
    19461944                                     drho_air_zw(k-1)               *          &
    19471945                                     ( zu(k+1) - zu(k-1)    )       /          &
    1948                                      ( km_sfc  + 1.0E-20_wp )   
     1946                                     ( km(k,j,i)  + 1.0E-20_wp )   
    19491947
    19501948          IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
     
    19651963          j = surf_def_h(1)%j(m)
    19661964          k = surf_def_h(1)%k(m)
    1967 
     1965!
     1966!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
     1967!--       and km are not on the same grid. Actually, a further
     1968!--       interpolation of km onto the u/v-grid is necessary. However, the
     1969!--       effect of this error is negligible.
    19681970          surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *          &
    19691971                                     drho_air_zw(k-1) *                        &
     
    19891991       DO  m = 1, surf_lsm_h%ns
    19901992
    1991           i = surf_lsm_h%i(m)           
     1993          i = surf_lsm_h%i(m)
    19921994          j = surf_lsm_h%j(m)
    19931995          k = surf_lsm_h%k(m)
    19941996!
    1995 !--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
     1997!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    19961998!--       and km are not on the same grid. Actually, a further
    19971999!--       interpolation of km onto the u/v-grid is necessary. However, the
    1998 !--       effect of this error is negligible.
    1999 !--       effect of this error is negligible.
    2000           km_sfc = surf_lsm_h%us(m) * kappa * surf_lsm_h%z_mo(m) /             &
    2001                    phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
    2002 
     2000!--       effect of this error is negligible.
    20032001          surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
    20042002                                        drho_air_zw(k-1)         *             &
    2005                                         ( zu(k+1) - zu(k-1)    ) /             &
    2006                                         ( km_sfc + 1.0E-20_wp ) 
     2003                                        ( zu(k+1)   - zu(k-1)    ) /           &
     2004                                        ( km(k,j,i) + 1.0E-20_wp ) 
    20072005          surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
    20082006                                        drho_air_zw(k-1)         *             &
    20092007                                        ( zu(k+1) - zu(k-1)    ) /             &
    2010                                         ( km_sfc  + 1.0E-20_wp )
     2008                                        ( km(k,j,i)  + 1.0E-20_wp )
    20112009
    20122010          IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
     
    20242022       DO  m = 1, surf_usm_h%ns
    20252023
    2026           i = surf_usm_h%i(m)           
     2024          i = surf_usm_h%i(m)
    20272025          j = surf_usm_h%j(m)
    20282026          k = surf_usm_h%k(m)
    20292027!
    2030 !--       Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
     2028!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    20312029!--       and km are not on the same grid. Actually, a further
    20322030!--       interpolation of km onto the u/v-grid is necessary. However, the
    2033 !--       effect of this error is negligible.
    2034           km_sfc = surf_usm_h%us(m) * kappa * surf_usm_h%z_mo(m) /             &
    2035                    phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
    2036 
     2031!--       effect of this error is negligible.
    20372032          surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
    20382033                                        drho_air_zw(k-1)         *             &
    20392034                                        ( zu(k+1) - zu(k-1)    ) /             &
    2040                                         ( km_sfc  + 1.0E-20_wp )
     2035                                        ( km(k,j,i)  + 1.0E-20_wp )
    20412036          surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
    20422037                                        drho_air_zw(k-1)         *             &
    20432038                                        ( zu(k+1) - zu(k-1)    ) /             &
    2044                                         ( km_sfc  + 1.0E-20_wp )
     2039                                        ( km(k,j,i)  + 1.0E-20_wp )
    20452040
    20462041          IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
     
    20572052
    20582053 END SUBROUTINE production_e_init
    2059 
    2060 
    2061 !------------------------------------------------------------------------------!
    2062 ! Description:
    2063 ! ------------
    2064 !> Calculates stability function for momentum
    2065 !>
    2066 !> @author Hauke Wurps
    2067 !------------------------------------------------------------------------------!
    2068  FUNCTION phi_m( zeta )
    2069 
    2070     IMPLICIT NONE
    2071 
    2072     REAL(wp)            :: phi_m         !< Value of the function
    2073     REAL(wp)            :: zeta          !< Stability parameter z/L
    2074 
    2075     REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
    2076     REAL(wp), PARAMETER :: b = -0.25_wp  !< constant
    2077     REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
    2078 
    2079     IF ( zeta < 0.0_wp )  THEN
    2080        phi_m = ( 1.0_wp - a * zeta )**b
    2081     ELSE
    2082        phi_m = 1.0_wp + c * zeta
    2083     ENDIF
    2084 
    2085  END FUNCTION phi_m
    20862054
    20872055
     
    46264594 SUBROUTINE tcm_diffusivities( var, var_reference )
    46274595
     4596    USE control_parameters,                                                    &
     4597        ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
     4598
     4599    USE surface_layer_fluxes_mod,                                              &
     4600        ONLY:  phi_m
     4601
     4602    USE surface_mod,                                                           &
     4603        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
     4604                surf_usm_h, surf_usm_v
     4605
     4606    INTEGER(iwp) ::  i                   !< loop index
     4607    INTEGER(iwp) ::  j                   !< loop index
     4608    INTEGER(iwp) ::  k                   !< loop index
     4609    INTEGER(iwp) ::  m                   !< loop index
     4610    INTEGER(iwp) ::  n                   !< loop index
     4611
    46284612    REAL(wp) ::  var_reference       !< reference temperature
    46294613
     
    46334617    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    46344618#endif
    4635 !
    4636 !-- Call default diffusivities routine. This is always used to calculate kh.
    4637     CALL tcm_diffusivities_default( var, var_reference )
    4638 !
    4639 !-- Call dynamic subgrid model to calculate km.
    4640     IF ( les_dynamic )  THEN
    4641        CALL tcm_diffusivities_dynamic
    4642     ENDIF
    4643 
    4644  END SUBROUTINE tcm_diffusivities
    4645 
    4646 
    4647 !------------------------------------------------------------------------------!
    4648 ! Description:
    4649 ! ------------
    4650 !> Computation of the turbulent diffusion coefficients for momentum and heat
    4651 !> according to Prandtl-Kolmogorov.
    4652 !------------------------------------------------------------------------------!
    4653  SUBROUTINE tcm_diffusivities_default( var, var_reference )
    4654  
    4655 
    4656     USE control_parameters,                                                    &
    4657         ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
    4658 
    4659     USE grid_variables,                                                        &
    4660         ONLY:  dx, dy
    4661 
    4662     USE statistics,                                                            &
    4663         ONLY :  rmask, sums_l_l
    4664 
    4665     USE surface_mod,                                                           &
    4666         ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
    4667                 surf_usm_h, surf_usm_v
    4668 
    4669     IMPLICIT NONE
    4670 
    4671     INTEGER(iwp) ::  i                   !< loop index
    4672     INTEGER(iwp) ::  j                   !< loop index
    4673     INTEGER(iwp) ::  k                   !< loop index
    4674     INTEGER(iwp) ::  m                   !< loop index
    4675     INTEGER(iwp) ::  n                   !< loop index
    4676     INTEGER(iwp) ::  omp_get_thread_num  !< opemmp function to get thread number
    4677     INTEGER(iwp) ::  sr                  !< statistic region
    4678     INTEGER(iwp) ::  tn                  !< thread number
    4679 
    4680     REAL(wp)     ::  flag                !< topography flag
    4681     REAL(wp)     ::  l                   !< mixing length
    4682     REAL(wp)     ::  ll                  !< adjusted mixing length
    4683     REAL(wp)     ::  var_reference       !< reference temperature
    4684 
    4685 #if defined( __nopointer )
    4686     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
    4687 #else
    4688     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    4689 #endif
    4690 
    4691 !
    4692 !-- Default thread number in case of one thread
    4693     tn = 0
    4694 
    4695 !
    4696 !-- Initialization for calculation of the mixing length profile
    4697     sums_l_l = 0.0_wp
    4698 
    4699 !
    4700 !-- Compute the turbulent diffusion coefficient for momentum
    4701     !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
    4702 !$  tn = omp_get_thread_num()
    47034619
    47044620!
    47054621!-- Introduce an optional minimum tke
    47064622    IF ( e_min > 0.0_wp )  THEN
    4707        !$OMP DO
     4623       !$OMP PARALLEL DO PRIVATE(i,j,k)
    47084624       DO  i = nxlg, nxrg
    47094625          DO  j = nysg, nyng
     
    47164632    ENDIF
    47174633
    4718     IF ( les_dynamic .OR. les_mw )  THEN
    4719        !$OMP DO
    4720        DO  i = nxlg, nxrg
    4721           DO  j = nysg, nyng
    4722              DO  k = nzb+1, nzt
    4723 
    4724                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4725 
    4726 !
    4727 !--             Determine the mixing length for LES closure
    4728                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
    4729 !
    4730 !--             Compute diffusion coefficients for momentum and heat
    4731                 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    4732                 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
    4733 !
    4734 !--             Summation for averaged profile (cf. flow_statistics)
    4735                 DO  sr = 0, statistic_regions
    4736                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
    4737                                                              * flag
    4738                 ENDDO
    4739 
    4740              ENDDO
    4741           ENDDO
    4742        ENDDO
    4743 
    4744     ELSEIF ( rans_tke_l )  THEN
    4745 
    4746        !$OMP DO
    4747        DO  i = nxlg, nxrg
    4748           DO  j = nysg, nyng
    4749              DO  k = nzb+1, nzt
    4750 
    4751                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4752 !
    4753 !--             Mixing length for RANS mode with TKE-l closure
    4754                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
    4755 !
    4756 !--             Compute diffusion coefficients for momentum and heat
    4757                 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    4758                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
    4759 !
    4760 !--             Summation for averaged profile (cf. flow_statistics)
    4761                 DO  sr = 0, statistic_regions
    4762                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
    4763                                                              * flag
    4764                 ENDDO
    4765 
    4766              ENDDO
    4767           ENDDO
    4768        ENDDO
    4769 
    4770     ELSEIF ( rans_tke_e )  THEN
    4771 
    4772        !$OMP DO
    4773        DO  i = nxlg, nxrg
    4774           DO  j = nysg, nyng
    4775              DO  k = nzb+1, nzt
    4776 
    4777                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4778 !
    4779 !--             Compute diffusion coefficients for momentum and heat
    4780                 km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
    4781                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
    4782 !
    4783 !--             Summation for averaged profile of mixing length (cf. flow_statistics)
    4784                 DO  sr = 0, statistic_regions
    4785                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
    4786                       c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
    4787                       ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
    4788                 ENDDO
    4789 
    4790              ENDDO
    4791           ENDDO
    4792        ENDDO
    4793 
     4634!
     4635!-- Call default diffusivities routine. This is always used to calculate kh.
     4636    CALL tcm_diffusivities_default( var, var_reference )
     4637!
     4638!-- Call dynamic subgrid model to calculate km.
     4639    IF ( les_dynamic )  THEN
     4640       CALL tcm_diffusivities_dynamic
    47944641    ENDIF
    47954642
    4796     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
    4797                                                 ! data output
    4798 !$OMP END PARALLEL
    4799 
    4800 !
    4801 !-- Set vertical boundary values (Neumann conditions both at upward- and
    4802 !-- downward facing walls. To set wall-boundary values, the surface data type
    4803 !-- is applied.
    4804 !-- Horizontal boundary conditions at vertical walls are not set because
    4805 !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary
    4806 !-- values of the diffusivities are not needed.
    4807     IF ( .NOT. rans_tke_e )  THEN
    4808 !
    4809 !--    Upward-facing
    4810        !$OMP PARALLEL DO PRIVATE( i, j, k )
    4811        DO  m = 1, bc_h(0)%ns
    4812           i = bc_h(0)%i(m)           
    4813           j = bc_h(0)%j(m)
    4814           k = bc_h(0)%k(m)
    4815           km(k-1,j,i) = km(k,j,i)
    4816           kh(k-1,j,i) = kh(k,j,i)
    4817        ENDDO
    4818 !
    4819 !--    Downward facing surfaces
    4820        !$OMP PARALLEL DO PRIVATE( i, j, k )
    4821        DO  m = 1, bc_h(1)%ns
    4822           i = bc_h(1)%i(m)           
    4823           j = bc_h(1)%j(m)
    4824           k = bc_h(1)%k(m)
    4825           km(k+1,j,i) = km(k,j,i)
    4826           kh(k+1,j,i) = kh(k,j,i)
    4827        ENDDO
    4828     ELSE
    4829 !
    4830 !--    Up- and downward facing surfaces
     4643!
     4644!-- Set boundary values (Neumann conditions at downward facing surfaces,
     4645!-- according to MOST at upward facing surfaces and vertical surfaces).
     4646!
     4647!-- Downward facing surfaces
     4648    !$OMP PARALLEL DO PRIVATE(i,j,k)
     4649    DO  m = 1, bc_h(1)%ns
     4650       i = bc_h(1)%i(m)
     4651       j = bc_h(1)%j(m)
     4652       k = bc_h(1)%k(m)
     4653       km(k+1,j,i) = km(k,j,i)
     4654       kh(k+1,j,i) = kh(k,j,i)
     4655    ENDDO
     4656
     4657!-- Upward facing surfaces
     4658!-- Default surfaces
     4659    n = 0
     4660    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4661    DO  m = 1, surf_def_h(n)%ns
     4662       i = surf_def_h(n)%i(m)
     4663       j = surf_def_h(n)%j(m)
     4664       k = surf_def_h(n)%k(m)
     4665       km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m) /       &
     4666                   phi_m( surf_def_h(n)%z_mo(m) / surf_def_h(n)%ol(m) )
     4667       kh(k,j,i) = 1.35_wp * km(k,j,i)
     4668    ENDDO
     4669!
     4670!-- Natural surfaces
     4671    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4672    DO  m = 1, surf_lsm_h%ns
     4673       i = surf_lsm_h%i(m)
     4674       j = surf_lsm_h%j(m)
     4675       k = surf_lsm_h%k(m)
     4676       km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /             &
     4677                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
     4678       kh(k,j,i) = 1.35_wp * km(k,j,i)
     4679    ENDDO
     4680!
     4681!-- Urban surfaces
     4682    !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4683    DO  m = 1, surf_usm_h%ns
     4684       i = surf_usm_h%i(m)
     4685       j = surf_usm_h%j(m)
     4686       k = surf_usm_h%k(m)
     4687       km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /             &
     4688                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
     4689       kh(k,j,i) = 1.35_wp * km(k,j,i)
     4690    ENDDO
     4691
     4692!
     4693!-- North-, south-, west and eastward facing surfaces
     4694    DO  n = 0, 3
     4695!
    48314696!--    Default surfaces
    4832        DO  n = 0, 1
    4833           DO  m = 1, surf_def_h(n)%ns
    4834              i = surf_def_h(n)%i(m)
    4835              j = surf_def_h(n)%j(m)
    4836              k = surf_def_h(n)%k(m)
    4837              km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m)
    4838              kh(k,j,i) = 1.35_wp * km(k,j,i)
    4839           ENDDO
    4840        ENDDO
    4841 !
    4842 !--    Upward facing surfaces
    4843 !--    Natural surfaces
    4844        DO  m = 1, surf_lsm_h%ns
    4845           i = surf_lsm_h%i(m)
    4846           j = surf_lsm_h%j(m)
    4847           k = surf_lsm_h%k(m)
    4848           km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m)
     4697       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4698       DO  m = 1, surf_def_v(n)%ns
     4699          i = surf_def_v(n)%i(m)
     4700          j = surf_def_v(n)%j(m)
     4701          k = surf_def_v(n)%k(m)
     4702          km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m) /    &
     4703                      phi_m( surf_def_v(n)%z_mo(m) / surf_def_v(n)%ol(m) )
    48494704          kh(k,j,i) = 1.35_wp * km(k,j,i)
    48504705       ENDDO
    48514706!
    4852 !--    Urban surfaces
    4853        DO  m = 1, surf_usm_h%ns
    4854           i = surf_usm_h%i(m)
    4855           j = surf_usm_h%j(m)
    4856           k = surf_usm_h%k(m)
    4857           km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m)
     4707!--    Natural surfaces
     4708       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4709       DO  m = 1, surf_lsm_v(n)%ns
     4710          i = surf_lsm_v(n)%i(m)
     4711          j = surf_lsm_v(n)%j(m)
     4712          k = surf_lsm_v(n)%k(m)
     4713          km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m) /    &
     4714                      phi_m( surf_lsm_v(n)%z_mo(m) / surf_lsm_v(n)%ol(m) )
    48584715          kh(k,j,i) = 1.35_wp * km(k,j,i)
    48594716       ENDDO
    4860 
    4861 !
    4862 !--    North-, south-, west and eastward facing surfaces
    4863        DO  n = 0, 3
    4864 !
    4865 !--       Default surfaces
    4866           DO  m = 1, surf_def_v(n)%ns
    4867              i = surf_def_v(n)%i(m)
    4868              j = surf_def_v(n)%j(m)
    4869              k = surf_def_v(n)%k(m)
    4870              km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
    4871              kh(k,j,i) = 1.35_wp * km(k,j,i)
    4872           ENDDO
    4873 !
    4874 !--       Natural surfaces
    4875           DO  m = 1, surf_lsm_v(n)%ns
    4876              i = surf_lsm_v(n)%i(m)
    4877              j = surf_lsm_v(n)%j(m)
    4878              k = surf_lsm_v(n)%k(m)
    4879              km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
    4880              kh(k,j,i) = 1.35_wp * km(k,j,i)
    4881           ENDDO
    4882 !
    4883 !--       Urban surfaces
    4884           DO  m = 1, surf_usm_v(n)%ns
    4885              i = surf_usm_v(n)%i(m)
    4886              j = surf_usm_v(n)%j(m)
    4887              k = surf_usm_v(n)%k(m)
    4888              km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
    4889              kh(k,j,i) = 1.35_wp * km(k,j,i)
    4890           ENDDO
     4717!
     4718!--    Urban surfaces
     4719       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
     4720       DO  m = 1, surf_usm_v(n)%ns
     4721          i = surf_usm_v(n)%i(m)
     4722          j = surf_usm_v(n)%j(m)
     4723          k = surf_usm_v(n)%k(m)
     4724          km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m) /    &
     4725                      phi_m( surf_usm_v(n)%z_mo(m) / surf_usm_v(n)%ol(m) )
     4726          kh(k,j,i) = 1.35_wp * km(k,j,i)
    48914727       ENDDO
    4892 
    4893        CALL exchange_horiz( km, nbgp )
    4894        CALL exchange_horiz( kh, nbgp )
    4895 
    4896     ENDIF
     4728    ENDDO
     4729
     4730    CALL exchange_horiz( km, nbgp )
     4731    CALL exchange_horiz( kh, nbgp )
    48974732
    48984733!
     
    49264761    ENDIF
    49274762
     4763 END SUBROUTINE tcm_diffusivities
     4764
     4765
     4766!------------------------------------------------------------------------------!
     4767! Description:
     4768! ------------
     4769!> Computation of the turbulent diffusion coefficients for momentum and heat
     4770!> according to Prandtl-Kolmogorov.
     4771!------------------------------------------------------------------------------!
     4772 SUBROUTINE tcm_diffusivities_default( var, var_reference )
     4773 
     4774    USE statistics,                                                            &
     4775        ONLY :  rmask, sums_l_l
     4776
     4777    IMPLICIT NONE
     4778
     4779    INTEGER(iwp) ::  i                   !< loop index
     4780    INTEGER(iwp) ::  j                   !< loop index
     4781    INTEGER(iwp) ::  k                   !< loop index
     4782    INTEGER(iwp) ::  omp_get_thread_num  !< opemmp function to get thread number
     4783    INTEGER(iwp) ::  sr                  !< statistic region
     4784    INTEGER(iwp) ::  tn                  !< thread number
     4785
     4786    REAL(wp)     ::  flag                !< topography flag
     4787    REAL(wp)     ::  l                   !< mixing length
     4788    REAL(wp)     ::  ll                  !< adjusted mixing length
     4789    REAL(wp)     ::  var_reference       !< reference temperature
     4790
     4791#if defined( __nopointer )
     4792    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
     4793#else
     4794    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
     4795#endif
     4796
     4797!
     4798!-- Default thread number in case of one thread
     4799    tn = 0
     4800
     4801!
     4802!-- Initialization for calculation of the mixing length profile
     4803    sums_l_l = 0.0_wp
     4804
     4805!
     4806!-- Compute the turbulent diffusion coefficient for momentum
     4807    !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
     4808!$  tn = omp_get_thread_num()
     4809
     4810    IF ( les_dynamic .OR. les_mw )  THEN
     4811       !$OMP DO
     4812       DO  i = nxlg, nxrg
     4813          DO  j = nysg, nyng
     4814             DO  k = nzb+1, nzt
     4815
     4816                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4817
     4818!
     4819!--             Determine the mixing length for LES closure
     4820                CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
     4821!
     4822!--             Compute diffusion coefficients for momentum and heat
     4823                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
     4824                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
     4825!
     4826!--             Summation for averaged profile (cf. flow_statistics)
     4827                DO  sr = 0, statistic_regions
     4828                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
     4829                                                             * flag
     4830                ENDDO
     4831
     4832             ENDDO
     4833          ENDDO
     4834       ENDDO
     4835
     4836    ELSEIF ( rans_tke_l )  THEN
     4837
     4838       !$OMP DO
     4839       DO  i = nxlg, nxrg
     4840          DO  j = nysg, nyng
     4841             DO  k = nzb+1, nzt
     4842
     4843                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4844!
     4845!--             Mixing length for RANS mode with TKE-l closure
     4846                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
     4847!
     4848!--             Compute diffusion coefficients for momentum and heat
     4849                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
     4850                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
     4851!
     4852!--             Summation for averaged profile (cf. flow_statistics)
     4853                DO  sr = 0, statistic_regions
     4854                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
     4855                                                             * flag
     4856                ENDDO
     4857
     4858             ENDDO
     4859          ENDDO
     4860       ENDDO
     4861
     4862    ELSEIF ( rans_tke_e )  THEN
     4863
     4864       !$OMP DO
     4865       DO  i = nxlg, nxrg
     4866          DO  j = nysg, nyng
     4867             DO  k = nzb+1, nzt
     4868
     4869                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4870!
     4871!--             Compute diffusion coefficients for momentum and heat
     4872                km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
     4873                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
     4874!
     4875!--             Summation for averaged profile of mixing length (cf. flow_statistics)
     4876                DO  sr = 0, statistic_regions
     4877                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
     4878                      c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
     4879                      ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
     4880                ENDDO
     4881
     4882             ENDDO
     4883          ENDDO
     4884       ENDDO
     4885
     4886    ENDIF
     4887
     4888    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
     4889                                                ! data output
     4890!$OMP END PARALLEL
     4891
    49284892 END SUBROUTINE tcm_diffusivities_default
    49294893
     
    49494913    USE arrays_3d,                                                             &
    49504914        ONLY:  ddzw, dzw, dd2zu, w, ug, vg
    4951 
    4952     USE control_parameters,                                                    &
    4953         ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
    49544915   
    49554916    USE grid_variables,                                                        &
    49564917        ONLY : ddx, ddy, dx, dy
    4957 
    4958     USE surface_mod,                                                           &
    4959         ONLY :  bc_h
    49604918
    49614919    IMPLICIT NONE
     
    49664924    INTEGER(iwp) ::  l           !< running index
    49674925    INTEGER(iwp) ::  m           !< running index
    4968     INTEGER(iwp) ::  n           !< running index
    49694926
    49704927    REAL(wp)     ::  dudx        !< Gradient of u-component in x-direction
     
    50374994    REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp))  !< constant
    50384995
    5039 !
    5040 !-- Introduce an optional minimum tke
    5041     IF ( e_min > 0.0_wp )  THEN
    5042        !$OMP DO
    5043        DO  i = nxlg, nxrg
    5044           DO  j = nysg, nyng
    5045              DO  k = nzb+1, nzt
    5046                 e(k,j,i) = MAX( e(k,j,i), e_min ) *                            &
    5047                         MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    5048              ENDDO
    5049           ENDDO
    5050        ENDDO
    5051     ENDIF
    50524996!
    50534997!-- velocities on grid centers:
     
    51815125       ENDDO
    51825126    ENDDO
    5183 !
    5184 !-- Set vertical boundary values (Neumann conditions both at upward- and
    5185 !-- downward facing walls. To set wall-boundary values, the surface data type
    5186 !-- is applied.
    5187 !-- Horizontal boundary conditions at vertical walls are not set because
    5188 !-- so far vertical surfaces require usage of a constant-flux layer where the
    5189 !-- boundary values of the diffusivities are not needed.
    5190 
    5191 !
    5192 !-- Upward-facing
    5193     !$OMP PARALLEL DO PRIVATE( i, j, k )
    5194     DO  m = 1, bc_h(0)%ns
    5195         i = bc_h(0)%i(m)           
    5196         j = bc_h(0)%j(m)
    5197         k = bc_h(0)%k(m)
    5198         km(k-1,j,i) = km(k,j,i)
    5199     ENDDO
    5200 !
    5201 !-- Downward facing surfaces
    5202     !$OMP PARALLEL DO PRIVATE( i, j, k )
    5203     DO  m = 1, bc_h(1)%ns
    5204         i = bc_h(1)%i(m)           
    5205         j = bc_h(1)%j(m)
    5206         k = bc_h(1)%k(m)
    5207         km(k+1,j,i) = km(k,j,i)
    5208     ENDDO
    5209 !
    5210 !-- Model top
    5211     !$OMP PARALLEL DO
    5212     DO  i = nxlg, nxrg
    5213        DO  j = nysg, nyng
    5214           km(nzt+1,j,i) = km(nzt,j,i)
    5215        ENDDO
    5216     ENDDO
    5217 !
    5218 !-- Set Neumann boundary conditions at the outflow boundaries in case of
    5219 !-- non-cyclic lateral boundaries
    5220     IF ( outflow_l )  THEN
    5221        km(:,:,nxl-1) = km(:,:,nxl)
    5222     ENDIF
    5223     IF ( outflow_r )  THEN
    5224        km(:,:,nxr+1) = km(:,:,nxr)
    5225     ENDIF
    5226     IF ( outflow_s )  THEN
    5227        km(:,nys-1,:) = km(:,nys,:)
    5228     ENDIF
    5229     IF ( outflow_n )  THEN
    5230        km(:,nyn+1,:) = km(:,nyn,:)
    5231     ENDIF
    52325127
    52335128 END SUBROUTINE tcm_diffusivities_dynamic
Note: See TracChangeset for help on using the changeset viewer.