source: palm/trunk/SOURCE/diffusivities.f90 @ 1

 SUBROUTINE diffusivities( theta )

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Log: diffusivities.f90,v $
! Revision 1.24  2006/04/26 12:16:26  raasch
! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
!
! Revision 1.23  2006/02/23 12:13:29  raasch
! nzb_2d replaced by nzb_s_inner/outer, option of a minimum tke
! quasi boundary conditions is now set for sums_l_l instead of sums_l
!
! Revision 1.22  2005/03/26 20:16:26  raasch
! Eddy diffusivities are also calculated for the ghost points in order to spare
! the communication otherwise needed for the exchange of these points.
! Neumann boundary conditions set at the outflow boundaries in case of
! non-cyclic lateral boundaries
!
! Revision 1.21  2004/01/30 10:23:12  raasch
! Scalar lower k index nzb replaced by 2d-array nzb_2d
!
! Revision 1.20  2003/03/14 13:40:09  raasch
! Loop optimization (SQRT(e) solved in a seperate loop, joining of loops)
!
! Revision 1.19  2003/03/12 16:28:56  raasch
! phi_m is initialized in declaration statement in order to avoid run time
! errors due to compiler problems on ibm and nec
!
! Revision 1.18  2002/12/19 14:32:23  raasch
! Correction of mixing length term (l(k)/ll(k)). The condition kh=3*km in
! the unstable case is now also exactly met in the wall adjustment region.
! Factor 0.7 in wall adjustment part replaced by variable
! wall_adjustment_factor, which is set to 1.8 in modules.f90.
! Missing OMP DO statement added.
!
! Revision 1.17  2002/06/11 12:56:49  raasch
! OpenMP directives added,
! array l changed from allocatable to automatic, because it must be on the stack
! in case of using OpenMP
!
! Revision 1.16  2001/08/21  08:30:34  08:30:34  raasch (Siegfried Raasch)
! Wall adjustment of mixing length to 0.7 z can be switched off
! Lower boundary condition for km and kh changed to Neumann
! 
! Revision 1.15  2001/03/30 07:22:09  raasch
! Near surface mixing length is limited to 0.7*zu,
! if adjust_mixing_length=T, l is modified at all gridpoints (previously at
! k=1 only),
! Translation of remaining German identifiers (variables, subroutines, etc.)
!
! Revision 1.14  2001/01/22 06:31:36  raasch
! Module test_variables removed
!
! Revision 1.13  2001/01/02 17:27:41  raasch
! -dpt_dz_d, dpt_dz_u
!
! Revision 1.12  2000/12/20 10:35:47  letzel
! All comments translated into English.
!
! Revision 1.11  2000/04/18 08:10:42  schroeter
! Revision 1.4 wieder rueckgaengig gemacht, das Stabilitaets-
! kriterium basiert nun wieder auf zentralen Differenzen
!
! Revision 1.9  2000/04/13 13:42:14  schroeter
! je nach Initialisierungsmodus (trocken/feucht) fliesst in
! Mischungswegberechnung pt oder vpt ein, wird durch entsprechende
! Variablenuebergabe geregelt
!
! Revision 1.8  99/02/17  09:17:39  09:17:39  raasch (Siegfried Raasch)
! Kriterium fuer reduzierten Mischungsweg im stabil geschichteten Fall enger
! gefasst (Schichtung muss sowohl oberhalb als auch unterhalb des betrachteten
! Gitterpunkts stabil sein)
! 
! Revision 1.7  1998/09/22 17:18:41  raasch
! Testweise cm = 0.4 bei k=1, aber vorerst auskommentiert
!
! Revision 1.6  1998/08/05 06:52:44  raasch
! Mischungsweganpassung an Prandtl nur bei k=1
!
! Revision 1.2  1998/03/11 11:49:26  raasch
! Anpassung des Mischungsweges an den Prandtlschen Mischungsweg moeglich
!
! Revision 1.1  1997/09/19 07:41:10  raasch
! Initial revision
!
!
! Description:
! ------------
! Computation of the turbulent diffusion coefficients for momentum and heat 
! according to Prandtl-Kolmogorov
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE control_parameters
    USE grid_variables
    USE indices
    USE pegrid
    USE statistics

    IMPLICIT NONE

    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn

    REAL    ::  dpt_dz, l_stable, phi_m = 1.0

    REAL    ::  theta(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)

    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e


!
!-- Default thread number in case of one thread
    tn = 0

!
!-- Initialization for calculation of the mixing length profile
    sums_l_l = 0.0

!
!-- Compute the turbulent diffusion coefficient for momentum
    !$OMP PARALLEL PRIVATE (dpt_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
!$  tn = omp_get_thread_num()

    !$OMP DO
    DO  i = nxl-1, nxr+1
       DO  j = nys-1, nyn+1

!
!--       Compute the Phi-function for a possible adaption of the mixing length
!--       to the Prandtl mixing length
          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
             IF ( rif(j,i) >= 0.0 )  THEN
                phi_m = 1.0 + 5.0 * rif(j,i)
             ELSE
                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
             ENDIF
          ENDIF
          
!
!--       Introduce an optional minimum tke
          IF ( e_min > 0.0 )  THEN
             DO  k = nzb_s_inner(j,i)+1, nzt
                e(k,j,i) = MAX( e(k,j,i), e_min )
             ENDDO
          ENDIF

!
!--       Calculate square root of e in a seperate loop, because it is used
!--       twice in the next loop (better vectorization)
          DO  k = nzb_s_inner(j,i)+1, nzt
             sqrt_e(k) = SQRT( e(k,j,i) )
          ENDDO

!
!--       Determine the mixing length
          DO  k = nzb_s_inner(j,i)+1, nzt
             dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
             IF ( dpt_dz > 0.0 ) THEN
                l_stable = 0.76 * sqrt_e(k) / &
                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
             ELSE
                l_stable = l_grid(k)
             ENDIF
!
!--          Adjustment of the mixing length
             IF ( wall_adjustment )  THEN
                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
             ELSE
                l(k)  = MIN( l_grid(k), l_stable )
                ll(k) = l_grid(k)
             ENDIF
             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
                l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
                ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
             ENDIF

!
!--          Compute diffusion coefficients for momentum and heat
             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)

          ENDDO

!
!--       Summation for averaged profile (cf. flow_statistics)
          DO  sr = 0, statistic_regions
             IF ( rmask(j,i,sr) /= 0.0 )  THEN
                DO  k = nzb_s_outer(j,i)+1, nzt
                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
                ENDDO
             ENDIF
          ENDDO

       ENDDO
    ENDDO

    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 bottom and top).
!-- Horizontal boundary conditions at vertical walls are not set because
!-- so far vertical walls require usage of a Prandtl-layer where the boundary
!-- values of the diffusivities are not needed
    !$OMP PARALLEL DO
    DO  i = nxl-1, nxr+1
       DO  j = nys-1, nyn+1
          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
          km(nzt+1,j,i)            = km(nzt,j,i)
          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,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 diffusivities