Changeset 51 for palm/trunk/SOURCE/diffusion_u.f90

Timestamp:

Mar 7, 2007 8:38:00 AM (18 years ago)

Author:

raasch

Message:

preliminary version, several changes to be explained later

File:

• palm/trunk/SOURCE/diffusion_u.f90 (modified) (8 diffs)

palm/trunk/SOURCE/diffusion_u.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
-! 
+! wall functions now include diabatic conditions, call of routine wall_fluxes
 ! 
 ! Former revisions:
@@ -30 +30 @@
 ! ------------
 ! Diffusion term of the u-component
+! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
+!        and slows down the speed on NEC about 5-10%
 !------------------------------------------------------------------------------!
 
@@ -55 +57 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
+       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+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(nzb:nzt+1)      ::  usvs
        REAL, DIMENSION(:,:),   POINTER ::  usws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -100 +103 @@
 !--          Wall functions at the north and south walls, respectively
              IF ( wall_u(j,i) /= 0.0 )  THEN
+
+!
+!--             Calculate the horizontal momentum flux u'v'
+                CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
+                                  usvs, 1.0, 0.0, 0.0, 0.0 )
+
                 DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-                   usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
-                   usvs   = -usvs * ABS( usvs )
                    kmyp_x = 0.25 * &
                             ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
@@ -132 +139 @@
                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
                                                   )                            &
-                                     + wall_u(j,i) * usvs                      &
+                                     + wall_u(j,i) * usvs(k)                   &
                                    ) * ddy
                 ENDDO
@@ -204 +211 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
+       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+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(nzb:nzt+1)      ::  usvs
        REAL, DIMENSION(:,:),   POINTER ::  usws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -246 +254 @@
 !--    Wall functions at the north and south walls, respectively
        IF ( wall_u(j,i) .NE. 0.0 )  THEN
+
+!
+!--       Calculate the horizontal momentum flux u'v'
+          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
+                            usvs, 1.0, 0.0, 0.0, 0.0 )
+
           DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-             usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
-             usvs   = -usvs * ABS( usvs )
              kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
              kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
@@ -276 +288 @@
                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
                                                   )                            &
-                                     + wall_u(j,i) * usvs                      &
+                                     + wall_u(j,i) * usvs(k)                   &
                                    ) * ddy
           ENDDO