Changeset 3121 for palm/trunk/SOURCE

Timestamp:

Jul 11, 2018 6:46:49 PM (6 years ago)

Author:

gronemeier

Message:

consider wall function for Km within production_e_init

File:

• palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (6 diffs)

palm/trunk/SOURCE/turbulence_closure_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! - created the function phi_m
+! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
+!   surfaces
+! 
+! 3120 2018-07-11 18:30:57Z gronemeier
 ! - changed tcm_diffusivities to tcm_diffusivities_default
 ! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
@@ -1977 +1982 @@
     INTEGER(iwp) ::  m      !< running index surface elements
 
+    REAL(wp)     :: km_sfc  !< km according to MOST
+
     IF ( constant_flux_layer )  THEN
 !
@@ -2001 +2008 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible.
+          km_sfc = surf_def_h(0)%us(m) * kappa * surf_def_h(0)%z_mo(m) /       &
+                   phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
+
           surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
                                      drho_air_zw(k-1)               *          &
                                      ( zu(k+1) - zu(k-1)    )       /          &
-                                     ( km(k,j,i)  + 1.0E-20_wp )
+                                     ( km_sfc  + 1.0E-20_wp )
           surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
                                      drho_air_zw(k-1)               *          &
                                      ( zu(k+1) - zu(k-1)    )       /          &
-                                     ( km(k,j,i)  + 1.0E-20_wp )   
+                                     ( km_sfc  + 1.0E-20_wp )   
 
           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
@@ -2060 +2070 @@
 !--       effect of this error is negligible. 
 !--       effect of this error is negligible. 
+          km_sfc = surf_lsm_h%us(m) * kappa * surf_lsm_h%z_mo(m) /             &
+                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
+
           surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km(k,j,i)  + 1.0E-20_wp )  
+                                        ( km_sfc  + 1.0E-20_wp )  
           surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
@@ -2091 +2104 @@
 !--       interpolation of km onto the u/v-grid is necessary. However, the
 !--       effect of this error is negligible. 
+          km_sfc = surf_usm_h%us(m) * kappa * surf_usm_h%z_mo(m) /             &
+                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
+
           surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 1.0E-20_wp )
           surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
                                         drho_air_zw(k-1)         *             &
                                         ( zu(k+1) - zu(k-1)    ) /             &
-                                        ( km(k,j,i)  + 1.0E-20_wp )
+                                        ( km_sfc  + 1.0E-20_wp )
 
           IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
@@ -2113 +2129 @@
 
  END SUBROUTINE production_e_init
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculates stability function for momentum
+!>
+!> @author Hauke Wurps
+!------------------------------------------------------------------------------!
+ FUNCTION phi_m( zeta ) 
+
+    IMPLICIT NONE 
+
+    REAL(wp)            :: phi_m         !< Value of the function
+    REAL(wp)            :: zeta          !< Stability parameter z/L
+
+    REAL(wp), PARAMETER :: a = 16.0_wp   !< constant
+    REAL(wp), PARAMETER :: b = -0.25_wp  !< constant
+    REAL(wp), PARAMETER :: c = 5.0_wp    !< constant
+
+    IF ( zeta < 0.0_wp )  THEN
+       phi_m = ( 1.0_wp - a * zeta )**b
+    ELSE
+       phi_m = 1.0_wp + c * zeta
+    ENDIF
+
+ END FUNCTION phi_m