Changeset 4476 for palm


Ignore:
Timestamp:
Mar 27, 2020 12:56:41 PM (15 months ago)
Author:
maronga
Message:

added a modified Deardorff subgrid-scale model

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r4473 r4476  
    2525! -----------------
    2626! $Id$
     27! Renamed variables for subgrids-scale model steering
     28!
     29! 4473 2020-03-25 21:04:07Z gronemeier
    2730! moved wall_adjustment_factor to turbulence_closure_mod
    2831!
     
    548551    CHARACTER (LEN=20)   ::  reference_state = 'initial_profile'          !< namelist parameter
    549552    CHARACTER (LEN=20)   ::  timestep_scheme = 'runge-kutta-3'            !< namelist parameter
    550     CHARACTER (LEN=20)   ::  turbulence_closure = 'Moeng_Wyngaard'        !< namelist parameter
     553    CHARACTER (LEN=20)   ::  turbulence_closure = '1.5-order'             !< namelist parameter
    551554    CHARACTER (LEN=40)   ::  topography = 'flat'                          !< namelist parameter
    552555    CHARACTER (LEN=64)   ::  host = '????'                                !< configuration identifier as given by palmrun option -c, ENVPAR namelist parameter provided by palmrun
     
    750753    LOGICAL ::  large_scale_subsidence = .FALSE.                 !< namelist parameter
    751754    LOGICAL ::  land_surface = .FALSE.                           !< use land surface model?
     755    LOGICAL ::  les_dai = .FALSE.                                !< use Dai et al. turbulence closure (modified 1.5-order closure) for LES mode. Shall replace the default 1.5-order closure
    752756    LOGICAL ::  les_dynamic = .FALSE.                            !< use dynamic subgrid model as turbulence closure for LES mode
    753     LOGICAL ::  les_mw = .FALSE.                                 !< use Moeng-Wyngaard turbulence closure for LES mode
     757    LOGICAL ::  les_default = .FALSE.                            !< use 1.5-order default turbulence closure for LES mode
    754758    LOGICAL ::  lsf_exception = .FALSE.                          !< use of lsf with buildings (temporary)?
    755759    LOGICAL ::  lsf_surf = .TRUE.                                !< use surface forcing (large scale forcing)?
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r4473 r4476  
    2525! -----------------
    2626! $Id$
     27! - added new LES closure after Dai et al. (2020), which provides much better grid
     28!   convergence in stable boundary layer runs. The implementation is experimental
     29!   at the moment and should be used with special care. The new SGS closure can be
     30!   switched on via turbulence_closure = '1.5-order-dai'
     31! - variable ml_wall_adjusted renamed to delta as it represents a grid size and
     32!   not a mixing length (see Equations 14 and 18 in Maronga et al. 2015, GMD)
     33! - nameing of turbulence closures revised:
     34!        'Moeng_Wyngaard' to '1.5-order'
     35!        'TKE-l' to 'tke-l'
     36!        'TKE-e' to 'tke-e'
     37! - LOGICAL steering variable renamed:
     38!        les_mw to les_default
     39!
     40! 4473 2020-03-25 21:04:07Z gronemeier
    2741! - rename l-grid to gridsize-geometric-mean
    2842!          l-wall to ml-wall-adjusted
     
    140154               initializing_actions, intermediate_timestep_count,             &
    141155               intermediate_timestep_count_max, km_constant,                  &
    142                les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
     156               les_dai, les_dynamic, les_default,                             &
     157               ocean_mode, plant_canopy, prandtl_number,                      &
    143158               pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
    144159               timestep_scheme, turbulence_closure,                           &
     
    210225
    211226    REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  ml_blackadar             !< mixing length according to Blackadar
    212     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ml_wall_adjusted         !< mixing length adjusted according to near-by walls
     227    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  delta                    !< grid size, possibly limited by wall adjustment factor
     228    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  distance_to_wall         !< distance to the surface/wall
    213229
    214230!
     
    584600    SELECT CASE ( TRIM( turbulence_closure ) )
    585601
     602       CASE ( '1.5-order' )
     603          les_default = .TRUE.
     604
     605       CASE ( '1.5-order-dai' )
     606          les_dai = .TRUE.
     607
    586608       CASE ( 'dynamic' )
    587609          les_dynamic = .TRUE.
    588610
    589        CASE ( 'Moeng_Wyngaard' )
    590           les_mw = .TRUE.
    591 
    592        CASE ( 'TKE-l' )
     611       CASE ( 'tke-l' )
    593612          rans_tke_l = .TRUE.
    594613          rans_mode = .TRUE.
    595614
    596        CASE ( 'TKE-e' )
     615       CASE ( 'tke-e' )
    597616          rans_tke_e = .TRUE.
    598617          rans_mode = .TRUE.
     
    11571176                   DO  j = nysg, nyng
    11581177                      DO  k = nzb+1, nzt
    1159                          km(k,j,i) = c_0 * SQRT( e_init ) * MIN( ml_wall_adjusted(k,j,i), &
     1178                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), &
    11601179                                                                 ml_blackadar(k) )
    11611180                      ENDDO
     
    12021221                DO  j = nysg, nyng
    12031222                   DO  k = nzb+1, nzt
    1204                       km(k,j,i) = c_0 * SQRT( e_init ) * ml_wall_adjusted(k,j,i)
     1223                      km(k,j,i) = c_0 * SQRT( e_init ) * delta(k,j,i)
    12051224                   ENDDO
    12061225                ENDDO
     
    14231442    REAL(wp) :: distance_edge_yz_down  !< distance of grid box center to its boundary along yz diagonal (down)
    14241443    REAL(wp) :: distance_edge_yz_up    !< distance of grid box center to its boundary along yz diagonal (up)
    1425     REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal (down)
     1444    REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal
    14261445    REAL(wp) :: distance_edge_xz_up    !< distance of grid box center to its boundary along xz diagonal (up)
    14271446    REAL(wp) :: distance_edge_xy       !< distance of grid box center to its boundary along xy diagonal
     
    14321451    REAL(wp), DIMENSION(nzb:nzt+1) ::  gridsize_geometric_mean  !< geometric mean of grid sizes dx, dy, dz
    14331452
    1434     ! ALLOCATE( gridsize_geometric_mean(nzb:nzt+1) )
    1435     ALLOCATE( ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1453
     1454    ALLOCATE( delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1455
     1456
     1457    IF ( les_dai )  THEN
     1458       ALLOCATE ( distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1459    ENDIF
    14361460!
    14371461!-- Initialize the mixing length in case of an LES-simulation
     
    14541478
    14551479!
    1456 !--    Initialize near-wall mixing length ml_wall_adjusted
     1480!--    Initialize the grid size variable (delta) and the distance to the wall, if needed
    14571481       DO  k = nzb, nzt+1
    1458           ml_wall_adjusted(k,:,:) = gridsize_geometric_mean(k)
     1482          delta(k,:,:) = gridsize_geometric_mean(k)
     1483
     1484!
     1485!--       If Dai et al. (2020) closure is used, the distance to the wall must be stored
     1486          IF ( les_dai )  THEN
     1487             distance_to_wall(k,:,:) = zu(k)
     1488          ENDIF
    14591489       ENDDO
    14601490
     
    15691599                         distance_corners_up = 9999999.9_wp
    15701600                      ENDIF
    1571 !
    1572 !--                   Adjust mixing length by wall-adjustment factor
    1573                       ml_wall_adjusted(k,j,i) = wall_adjustment_factor * MIN(  &
     1601
     1602!
     1603!--                   Calculate the minimum distance from the wall and store it
     1604!--                   temporarily in the array delta
     1605                      delta(k,j,i) = MIN(                                      &
    15741606                         distance_up, distance_down, distance_ns, distance_lr, &
    15751607                         distance_edge_yz_down, distance_edge_yz_up,           &
     
    15781610                         distance_corners_down, distance_corners_up )
    15791611
     1612!
     1613!--                   If Dai et al. (2020) closure is used, the distance to the wall
     1614!--                   must be permanently stored
     1615                      IF ( les_dai  .AND.  delta(k,j,i) /= 9999999.9_wp )  THEN
     1616                         distance_to_wall(k,j,i) = delta(k,j,i)
     1617                      ENDIF
     1618
     1619!
     1620!--                   Close to the surface, the maximum mixing length is limited to 1.8 z
     1621                      delta(k,j,i) = wall_adjustment_factor * delta(k,j,i)
     1622
    15801623                  ENDIF  !if grid point belongs to atmosphere
    15811624
    1582                   ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &
    1583                                                  gridsize_geometric_mean(k) )
     1625
     1626
     1627!
     1628!--               The grid size (delta) is defined as the the minimum of the distance to
     1629!--               the nearest wall * 1.8 and the geometric mean grid size.
     1630                  delta(k,j,i) = MIN( delta(k,j,i), gridsize_geometric_mean(k) )
     1631
    15841632
    15851633                ENDDO  !k loop
     
    16221670       ENDDO
    16231671
    1624        ml_wall_adjusted(nzb,:,:) = ml_blackadar(nzb)
    1625        ml_wall_adjusted(nzt+1,:,:) = ml_blackadar(nzt+1)
     1672       delta(nzb,:,:) = ml_blackadar(nzb)
     1673       delta(nzt+1,:,:) = ml_blackadar(nzt+1)
    16261674!
    16271675!--    Limit mixing length to either nearest wall or Blackadar mixing length.
     
    16341682          radius = ml_blackadar(k)  ! radius within walls are searched
    16351683!
    1636 !--       Set ml_wall_adjusted to its default maximum value (ml_blackadar)
    1637           ml_wall_adjusted(k,:,:) = radius
     1684!--       Set delta to its default maximum value (ml_blackadar)
     1685          delta(k,:,:) = radius
    16381686
    16391687!
     
    17301778!--                            Search for walls within octant (+++)
    17311779                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1732                                ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
     1780                               delta(k,j,i) = MIN( delta(k,j,i),             &
    17331781                                       shortest_distance( vic_yz, .TRUE., ii ) )
    17341782!
     
    17381786!--                            shortest_distance").
    17391787                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1740                                ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
     1788                               delta(k,j,i) = MIN( delta(k,j,i),             &
    17411789                                       shortest_distance( vic_yz, .TRUE., ii ) )
    17421790
     
    17451793!--                         Search for walls within octant (+--)
    17461794                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1747                             ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
     1795                            delta(k,j,i) = MIN( delta(k,j,i),                &
    17481796                                      shortest_distance( vic_yz, .FALSE., ii ) )
    17491797!
    17501798!--                         Search for walls within octant (++-)
    17511799                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1752                             ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
     1800                            delta(k,j,i) = MIN( delta(k,j,i),                &
    17531801                                      shortest_distance( vic_yz, .FALSE., ii ) )
    17541802!
    17551803!--                         Reduce search along x by already found distance
    1756                             dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
     1804                            dist_dx = CEILING( delta(k,j,i) / dx )
    17571805
    17581806                         ENDDO
     
    17671815!--                            Search for walls within octant (-++)
    17681816                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1769                                ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
     1817                               delta(k,j,i) = MIN( delta(k,j,i),             &
    17701818                                      shortest_distance( vic_yz, .TRUE., -ii ) )
    17711819!
     
    17751823!--                            shortest_distance").
    17761824                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1777                                ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
     1825                               delta(k,j,i) = MIN( delta(k,j,i),             &
    17781826                                      shortest_distance( vic_yz, .TRUE., -ii ) )
    17791827
     
    17821830!--                         Search for walls within octant (---)
    17831831                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1784                             ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
     1832                            delta(k,j,i) = MIN( delta(k,j,i),                &
    17851833                                     shortest_distance( vic_yz, .FALSE., -ii ) )
    17861834!
    17871835!--                         Search for walls within octant (-+-)
    17881836                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1789                             ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
     1837                            delta(k,j,i) = MIN( delta(k,j,i),                &
    17901838                                     shortest_distance( vic_yz, .FALSE., -ii ) )
    17911839!
    17921840!--                         Reduce search along x by already found distance
    1793                             dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
     1841                            dist_dx = CEILING( delta(k,j,i) / dx )
    17941842
    17951843                         ENDDO
     
    17991847                   ELSE  !Check if (i,j,k) belongs to atmosphere
    18001848
    1801                       ml_wall_adjusted(k,j,i) = ml_blackadar(k)
     1849                      delta(k,j,i) = ml_blackadar(k)
    18021850
    18031851                   ENDIF
     
    18181866
    18191867!
    1820 !-- Set lateral boundary conditions for ml_wall_adjusted
    1821     CALL exchange_horiz( ml_wall_adjusted, nbgp )
    1822 
    1823     !$ACC ENTER DATA COPYIN(ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     1868!-- Set lateral boundary conditions for delta
     1869    CALL exchange_horiz( delta, nbgp )
     1870    !$ACC ENTER DATA COPYIN(delta(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     1871
     1872    IF ( les_dai )  THEN
     1873       CALL exchange_horiz( distance_to_wall, nbgp )
     1874       !$ACC ENTER DATA COPYIN(distance_to_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     1875    ENDIF
     1876
     1877
    18241878
    18251879    CONTAINS
     
    40324086!
    40334087!-- Calculate the dissipation
    4034     IF ( les_dynamic .OR. les_mw )  THEN
     4088    IF( les_default  .OR.  les_dynamic )  THEN
    40354089
    40364090       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    40374091       !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
    40384092       !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
    4039        !$ACC PRESENT(dd2zu, ml_wall_adjusted)
     4093       !$ACC PRESENT(dd2zu, delta)
    40404094       DO  i = nxl, nxr
    40414095          DO  j = nys, nyn
     
    40534107                   ENDIF
    40544108                ELSE
    4055                    ml_stratification(k) = ml_wall_adjusted(k,j,i)
     4109                   ml_stratification(k) = delta(k,j,i)
    40564110                ENDIF
    40574111
    40584112             ENDDO
    40594113
     4114!
     4115!--          ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by
     4116!--          the Intel18 compiler. This behaviour might change for future compiler versions.
    40604117             !$ACC LOOP PRIVATE(k)
    40614118             !DIR$ IVDEP
    40624119             DO  k = nzb+1, nzt
    40634120
    4064                 ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
    4065 
    4066                 diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
     4121                ml = MIN( delta(k,j,i), ml_stratification(k) )
     4122
     4123                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )           &
    40674124                              * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
    40684125                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     
    40734130       ENDDO
    40744131
     4132    ELSEIF ( les_dai )  THEN
     4133
     4134       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4135       !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
     4136       !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
     4137       !$ACC PRESENT(dd2zu, delta)
     4138       DO  i = nxl, nxr
     4139          DO  j = nys, nyn
     4140             !$ACC LOOP PRIVATE(k)
     4141             DO  k = nzb+1, nzt
     4142                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4143                IF ( dvar_dz > 0.0_wp ) THEN
     4144!
     4145!--                The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
     4146!--                the stratification term. 1E-5 is added as l is zero at the beginning of
     4147!--                the simulation.
     4148                   IF ( use_single_reference_value )  THEN
     4149                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) &
     4150                                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4151                   ELSE
     4152                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) ) &
     4153                                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4154                   ENDIF
     4155                   ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
     4156                                                   + 1.0_wp / ml_stratification(k) )
     4157                ELSE
     4158                   ml_stratification(k) = delta(k,j,i)
     4159                ENDIF
     4160
     4161             ENDDO
     4162
     4163!
     4164!--          ATTENTION: Don't merge the following loop with the previous one, because this would prohibit proper vectorization by
     4165!--          the Intel18 compiler. This behaviour might change for future compiler versions.
     4166             !$ACC LOOP PRIVATE(k)
     4167             !DIR$ IVDEP
     4168             DO  k = nzb+1, nzt
     4169
     4170                ml          = ml_stratification(k)
     4171                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
     4172                              * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
     4173                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4174
     4175             ENDDO
     4176
     4177          ENDDO
     4178       ENDDO
    40754179    ELSEIF ( rans_tke_l )  THEN
    40764180
     
    40784182      !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) &
    40794183      !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
    4080       !$ACC PRESENT(dd2zu, ml_blackadar, ml_wall_adjusted)
     4184      !$ACC PRESENT(dd2zu, ml_blackadar, delta)
    40814185       DO  i = nxl, nxr
    40824186          DO  j = nys, nyn
     
    40844188!--          Calculate Richardson-flux number
    40854189             IF ( use_single_reference_value )  THEN
     4190
    40864191                !$ACC LOOP PRIVATE(k)
    40874192                DO  k = nzb+1, nzt
     
    40964201                ENDDO
    40974202             ELSE
     4203!
     4204!--             ATTENTION: Don't merge the following loops with the previous one, because this would prohibit proper vectorization
     4205!--             by the Intel18 compiler. This behaviour might change for future compiler versions.
    40984206                !$ACC LOOP PRIVATE(k)
    40994207                DO  k = nzb+1, nzt
     
    41244232
    41254233                diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
    4126                             / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
     4234                            / MIN( ml_stratification(k), delta(k,j,i) )                    &
    41274235                            * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    41284236
     
    42434351!-- Calculate the mixing length and dissipation...
    42444352!-- ...in case of LES
    4245     IF ( les_dynamic .OR. les_mw )  THEN
     4353    IF ( les_default  .OR.  les_dynamic )  THEN
    42464354
    42474355       DO  k = nzb+1, nzt
     
    42564364             ENDIF
    42574365          ELSE
    4258              ml_stratification(k) = ml_wall_adjusted(k,j,i)
     4366             ml_stratification(k) = delta(k,j,i)
    42594367          ENDIF
    42604368       ENDDO
     
    42624370       !DIR$ IVDEP
    42634371       DO  k = nzb+1, nzt
    4264           ml  = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
    4265 
    4266           diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
     4372          ml  = MIN( delta(k,j,i), ml_stratification(k) )
     4373
     4374          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
    42674375                        * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
    42684376                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    42694377       ENDDO
     4378
     4379    ELSEIF ( les_dai )  THEN
     4380
     4381       DO  k = nzb+1, nzt
     4382          dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4383          IF ( dvar_dz > 0.0_wp ) THEN
     4384!
     4385!--          The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
     4386!--          the stratification term. 1E-5 is added as l is zero at the beginning of
     4387!--          the simulation.
     4388             IF ( use_single_reference_value )  THEN
     4389                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
     4390                                       / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4391             ELSE
     4392                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
     4393                                       / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4394             ENDIF
     4395
     4396             ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
     4397                                             + 1.0_wp / ml_stratification(k) )
     4398
     4399          ELSE
     4400             ml_stratification(k) = delta(k,j,i)
     4401          ENDIF
     4402       ENDDO
     4403
     4404       !DIR$ IVDEP
     4405       DO  k = nzb+1, nzt
     4406          ml  = ml_stratification(k)
     4407
     4408          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / delta(k,j,i) )                      &
     4409                        * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
     4410                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4411       ENDDO
     4412
     4413
     4414
    42704415
    42714416!
     
    43104455       DO  k = nzb+1, nzt
    43114456          diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
    4312                       / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
     4457                      / MIN( ml_stratification(k), delta(k,j,i) )                    &
    43134458                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    43144459       ENDDO
     
    47304875!$  tn = omp_get_thread_num()
    47314876
    4732     IF ( les_dynamic .OR. les_mw )  THEN
     4877    IF ( les_default  .OR.  les_dynamic )  THEN
    47334878       !$OMP DO
    47344879       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    47354880       !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
    4736        !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, ml_wall_adjusted) &
     4881       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) &
    47374882       !$ACC PRESENT(kh, km, sums_l_l, rmask)
    47384883       DO  i = nxlg, nxrg
     
    47454890!--             due to errors when using OpenACC directives. The execution
    47464891!--             crashes reliably if a subroutine is called at this point (the
    4747 !--             reasong for this behaviour is unknown, however).
     4892!--             reason for this behaviour is unknown, however).
    47484893                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    47494894                IF ( dvar_dz > 0.0_wp ) THEN
     
    47564901                   ENDIF
    47574902                ELSE
    4758                    ml_stratification(k) = ml_wall_adjusted(k,j,i)
     4903                   ml_stratification(k) = delta(k,j,i)
    47594904                ENDIF
    47604905
     
    47654910             DO  k = nzb+1, nzt
    47664911
    4767                 ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )         &
     4912                ml = MIN( delta(k,j,i), ml_stratification(k) )                    &
    47684913                   * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    47694914                ml_local_profile(k) = ml
     
    47714916!--             Compute diffusion coefficients for momentum and heat
    47724917                km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) )
    4773                 kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / ml_wall_adjusted(k,j,i) ) * km(k,j,i)
     4918                kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / delta(k,j,i) ) * km(k,j,i)
    47744919
    47754920             ENDDO
     
    47864931       ENDDO
    47874932
     4933    ELSEIF ( les_dai )  THEN
     4934       !$OMP DO
     4935       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4936       !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
     4937       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, delta) &
     4938       !$ACC PRESENT(kh, km, sums_l_l, rmask)
     4939       DO  i = nxlg, nxrg
     4940          DO  j = nysg, nyng
     4941             !$ACC LOOP PRIVATE(k)
     4942             DO  k = nzb+1, nzt
     4943
     4944                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4945                IF ( dvar_dz > 0.0_wp ) THEN
     4946!
     4947!--                The mixing length is calculated as 1/l = 1/(kappa*z) + 1/Lb, where Lb is
     4948!--                the stratification term. 1E-5 is added as l is zero at the beginning of
     4949!--                the simulation.
     4950                   IF ( use_single_reference_value )  THEN
     4951                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
     4952                                             / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4953                   ELSE
     4954                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )                           &
     4955                                             / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4956                   ENDIF
     4957
     4958                   ml_stratification(k) = 1.0_wp / ( 1.0_wp / ( kappa * distance_to_wall(k,j,i) ) &
     4959                                        + 1.0_wp / ml_stratification(k) )
     4960                ELSE
     4961                   ml_stratification(k) = delta(k,j,i)
     4962                ENDIF
     4963
     4964             ENDDO
     4965
     4966             !$ACC LOOP PRIVATE(k)
     4967             !DIR$ IVDEP
     4968             DO  k = nzb+1, nzt
     4969
     4970                ml_local_profile(k) = ml_stratification(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4971                ml  = ml_local_profile(k)
     4972!
     4973!--             Compute diffusion coefficients for momentum and heat
     4974                km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) )
     4975
     4976                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4977                IF ( dvar_dz > 0.0_wp ) THEN
     4978                   kh(k,j,i) = km(k,j,i)
     4979                ELSE
     4980                   kh(k,j,i) = 3.0_wp * km(k,j,i)
     4981                ENDIF
     4982
     4983             ENDDO
     4984!
     4985!--          Summation for averaged profile (cf. flow_statistics)
     4986             !$ACC LOOP PRIVATE(sr, k)
     4987             DO  sr = 0, statistic_regions
     4988                DO  k = nzb+1, nzt
     4989                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
     4990                ENDDO
     4991             ENDDO
     4992
     4993          ENDDO
     4994       ENDDO
     4995
    47884996    ELSEIF ( rans_tke_l )  THEN
    47894997
     
    47914999       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    47925000       !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, rif) &
    4793        !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, ml_wall_adjusted, ml_blackadar) &
     5001       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, delta, ml_blackadar) &
    47945002       !$ACC PRESENT(kh, km, sums_l_l, rmask)
    47955003       DO  i = nxlg, nxrg
     
    48395047             !DIR$ IVDEP
    48405048             DO  k = nzb+1, nzt
    4841                 ml_local_profile(k) = MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) ) &
     5049                ml_local_profile(k) = MIN( ml_stratification(k), delta(k,j,i) ) &
    48425050                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    48435051                km(k,j,i) = c_0 * ml_local_profile(k) * SQRT( e(k,j,i) )
     
    49935201    REAL(wp)     ::  ldsd        !< sum: ld_ij*sd_ij
    49945202
    4995     REAL(wp)     ::  delta       !< grid size
     5203    REAL(wp)     ::  delta_max   !< maximum of the grid spacings
    49965204    REAL(wp)     ::  cst         !< c*
    49975205    REAL(wp)     ::  cstnust_t   !< product c*nu*
     
    51195327!--          The model was only tested for an isotropic grid. The following
    51205328!--          expression was a recommendation of Stefan Heinz.
    5121              delta = MAX( dx, dy, dzw(k) )
     5329             delta_max = MAX( dx, dy, dzw(k) )
    51225330
    51235331             IF ( lnn <= 0.0_wp ) THEN
     
    51255333             ELSE
    51265334                cst = cstnust_t /                                              &
    5127                    ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
     5335                   ( 4.0_wp * delta_max * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
    51285336             ENDIF
    51295337
     
    51315339!--          Calculate border according to Mokhtarpoor and Heinz (2017)
    51325340             cst_max = fac_cmax * SQRT( e(k,j,i) ) /                           &
    5133                        ( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
     5341                       ( delta_max * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
    51345342
    51355343             IF ( ABS( cst ) > cst_max )  THEN
     
    51375345             ENDIF
    51385346
    5139              km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
     5347             km(k,j,i) = cst * delta_max * SQRT( e(k,j,i) ) * flag
    51405348
    51415349          ENDDO
Note: See TracChangeset for help on using the changeset viewer.