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

 MODULE diffusion_w_mod

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Log: diffusion_w.f90,v $
! Revision 1.12  2006/02/23 10:38:03  raasch
! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
! +z0 in argument list
! WARNING: loops containing the MAX function are still not properly vectorized!
!
! Revision 1.11  2005/03/26 20:11:05  raasch
! Additional damping layer at the outflow in case of non-cyclic lateral
! boundaries, additional arguments km_damp_x, km_damp_y
!
! Revision 1.10  2004/01/30 10:22:32  raasch
! Scalar lower k index nzb replaced by 2d-array nzb_2d
!
! Revision 1.9  2003/03/12 16:26:49  raasch
! Full code replaced in the call for all gridpoints instead of calling the
! _ij version (required by NEC, because otherwise no vectorization)
!
! Revision 1.8  2002/06/11 12:54:16  raasch
! Former subroutine changed to a module which allows to be called for all grid
! points of a single vertical column with index i,j or for all grid points by
! using function overloading.
!
! Revision 1.7  2001/03/30 07:12:35  raasch
! Translation of remaining German identifiers (variables, subroutines, etc.)
!
! Revision 1.6  2001/01/22 06:20:12  raasch
! Module test_variables removed
!
! Revision 1.5  2000/07/03 12:58:40  raasch
! dummy arguments, whose corresponding actual arguments are pointers,
! are now also defined as pointers,
! all comments translated into English
!
! Revision 1.4  1998/07/06 12:12:24  raasch
! + USE test_variables
!
! Revision 1.3  1997/09/12 07:24:58  raasch
! Leerzeilen mussten entfernt werden
!
! Revision 1.2  1997/09/12 06:44:00  raasch
! HP-Compiler erfordert & am Beginn von Fortsetzungszeilen
!
! Revision 1.1  1997/09/12 06:24:11  raasch
! Initial revision
!
!
! Description:
! ------------
! Diffusion term of the w-component
!------------------------------------------------------------------------------!

    PRIVATE
    PUBLIC diffusion_w

    INTERFACE diffusion_w
       MODULE PROCEDURE diffusion_w
       MODULE PROCEDURE diffusion_w_ij
    END INTERFACE diffusion_w

 CONTAINS


!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
                            w, z0 )

       USE control_parameters
       USE grid_variables
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
                   kmyp_z, wsus, wsvs
       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
                   km_damp_y(nys-1:nyn+1)
       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w


       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_w_outer(j,i)+1, nzt-1
!
!--             Interpolate eddy diffusivities on staggered gridpoints
                kmxp_x = 0.25 * &
                         ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
                kmxm_x = 0.25 * &
                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
                kmxp_z = kmxp_x
                kmxm_z = kmxm_x
                kmyp_y = 0.25 * &
                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
                kmym_y = 0.25 * &
                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
                kmyp_z = kmyp_y
                kmym_z = kmym_y
!
!--             Increase diffusion at the outflow boundary in case of
!--             non-cyclic lateral boundaries. Damping is only needed for
!--             velocity components parallel to the outflow boundary in
!--             the direction normal to the outflow boundary.
                IF ( bc_lr /= 'cyclic' )  THEN
                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
                ENDIF
                IF ( bc_ns /= 'cyclic' )  THEN
                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
                   kmym_y = MAX( kmym_y, km_damp_y(j) )
                ENDIF

                tend(k,j,i) = tend(k,j,i)                                      &
                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
                      &   ) * ddx                                              &
                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
                      &   ) * ddy                                              &
                      & + 2.0 * (                                              &
                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
                      &         ) * ddzu(k+1)
             ENDDO

!
!--          Wall functions at all vertical walls, where necessary
             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
                   IF ( wall_w_x(j,i) /= 0.0 )  THEN
                      wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
                      wsus = -wsus * ABS( wsus )
                   ELSE
                      wsus = 0.0
                   ENDIF
                   IF ( wall_w_y(j,i) /= 0.0 )  THEN
                      wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
                      wsvs = -wsvs * ABS( wsvs )
                   ELSE
                      wsvs = 0.0
                   ENDIF
!
!--                Interpolate eddy diffusivities on staggered gridpoints
                   kmxp_x = 0.25 * &
                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
                   kmxm_x = 0.25 * &
                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
                   kmxp_z = kmxp_x
                   kmxm_z = kmxm_x
                   kmyp_y = 0.25 * &
                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
                   kmym_y = 0.25 * &
                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
                   kmyp_z = kmyp_y
                   kmym_z = kmym_y
!
!--                Increase diffusion at the outflow boundary in case of
!--                non-cyclic lateral boundaries. Damping is only needed for
!--                velocity components parallel to the outflow boundary in
!--                the direction normal to the outflow boundary.
                   IF ( bc_lr /= 'cyclic' )  THEN
                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
                   ENDIF
                   IF ( bc_ns /= 'cyclic' )  THEN
                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
                      kmym_y = MAX( kmym_y, km_damp_y(j) )
                   ENDIF

                   tend(k,j,i) = tend(k,j,i)                                   &
                                 + (   fwxp(j,i) * (                           &
                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
                                                   )                           &
                                     - fwxm(j,i) * (                           &
                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
                                                   )                           &
                                     + wall_w_x(j,i) * wsus                    &
                                   ) * ddx                                     &
                                 + (   fwyp(j,i) * (                           &
                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
                                                   )                           &
                                     - fwym(j,i) * (                           &
                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
                                                   )                           &
                                     + wall_w_y(j,i) * wsvs                    &
                                   ) * ddy                                     &
                                 + 2.0 * (                                     &
                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
                                         ) * ddzu(k+1)
                ENDDO
             ENDIF

          ENDDO
       ENDDO

    END SUBROUTINE diffusion_w


!------------------------------------------------------------------------------!
! Call for grid point i,j
!------------------------------------------------------------------------------!
    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
                               tend, u, v, w, z0 )

       USE control_parameters
       USE grid_variables
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
                   kmyp_z, wsus, wsvs
       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt), km_damp_x(nxl-1:nxr+1), &
                   km_damp_y(nys-1:nyn+1)
       REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w


       DO  k = nzb_w_outer(j,i)+1, nzt-1
!
!--       Interpolate eddy diffusivities on staggered gridpoints
          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
          kmxp_z = kmxp_x
          kmxm_z = kmxm_x
          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
          kmyp_z = kmyp_y
          kmym_z = kmym_y
!
!--       Increase diffusion at the outflow boundary in case of non-cyclic
!--       lateral boundaries. Damping is only needed for velocity components
!--       parallel to the outflow boundary in the direction normal to the
!--       outflow boundary.
          IF ( bc_lr /= 'cyclic' )  THEN
             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
          ENDIF
          IF ( bc_ns /= 'cyclic' )  THEN
             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
             kmym_y = MAX( kmym_y, km_damp_y(j) )
          ENDIF

          tend(k,j,i) = tend(k,j,i)                                            &
                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
                      &   ) * ddx                                              &
                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
                      &   ) * ddy                                              &
                      & + 2.0 * (                                              &
                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
                      &         ) * ddzu(k+1)
       ENDDO

!
!--    Wall functions at all vertical walls, where necessary
       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
             IF ( wall_w_x(j,i) /= 0.0 )  THEN
                wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
                wsus = -wsus * ABS( wsus )
             ELSE
                wsus = 0.0
             ENDIF
             IF ( wall_w_y(j,i) /= 0.0 )  THEN
                wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
                wsvs = -wsvs * ABS( wsvs )
             ELSE
                wsvs = 0.0
             ENDIF
!
!--          Interpolate eddy diffusivities on staggered gridpoints
             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
             kmxp_z = kmxp_x
             kmxm_z = kmxm_x
             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
             kmyp_z = kmyp_y
             kmym_z = kmym_y
!
!--          Increase diffusion at the outflow boundary in case of
!--          non-cyclic lateral boundaries. Damping is only needed for
!--          velocity components parallel to the outflow boundary in
!--          the direction normal to the outflow boundary.
             IF ( bc_lr /= 'cyclic' )  THEN
                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
             ENDIF
             IF ( bc_ns /= 'cyclic' )  THEN
                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
                kmym_y = MAX( kmym_y, km_damp_y(j) )
             ENDIF

             tend(k,j,i) = tend(k,j,i)                                         &
                                 + (   fwxp(j,i) * (                           &
                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
                                                   )                           &
                                     - fwxm(j,i) * (                           &
                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
                                                   )                           &
                                     + wall_w_x(j,i) * wsus                    &
                                   ) * ddx                                     &
                                 + (   fwyp(j,i) * (                           &
                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
                                                   )                           &
                                     - fwym(j,i) * (                           &
                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
                                                   )                           &
                                     + wall_w_y(j,i) * wsvs                    &
                                   ) * ddy                                     &
                                 + 2.0 * (                                     &
                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
                                         ) * ddzu(k+1)
          ENDDO
       ENDIF

    END SUBROUTINE diffusion_w_ij

 END MODULE diffusion_w_mod