minor code update: removing zeta from output parameters in neutral case

srcF/sfx_esm.f90

@@ -212,8 +212,9 @@ contains
         else if (Rib > -0.001) then
             ! --- nearly neutral [-0.001, 0] block
 
-            call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+            call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
 
+            zeta = 0.0
             phi_m = 1.0
             phi_h = 1.0 / Pr_t_0_inv
         else
@@ -268,19 +269,18 @@ contains
 
     ! universal functions
     ! --------------------------------------------------------------------------------
-    subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+    subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
         !> @brief universal functions (momentum) & (heat): neutral case
         ! ----------------------------------------------------------------------------
         real, intent(out) :: psi_m, psi_h   !> universal functions
-        real, intent(out) :: zeta           !> = z/L
 
         real, intent(in) :: h0_m, h0_t      !> = z/z0_m, z/z0_h
         real, intent(in) :: B               !> = log(z0_m / z0_h)
         ! ----------------------------------------------------------------------------
 
-        zeta = 0.0
         psi_m = log(h0_m)
         psi_h = log(h0_t) / Pr_t_0_inv
+        !*: this looks redundant z0_t = z0_m in case |B| ~ 0
         if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
 
     end subroutine

srcF/sfx_log.f90

@@ -99,7 +99,6 @@ contains
         real u_dyn0             !> dynamic velocity in neutral conditions [m/s]
         real Re                 !> roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
 
-        real zeta               !> = z/L [n/d]
         real Rib                !> bulk Richardson number
 
         real psi_m, psi_h       !> universal functions (momentum) & (heat) [n/d]
@@ -172,7 +171,7 @@ contains
 
         ! --- get the fluxes
         ! ----------------------------------------------------------------------------
-        call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+        call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
         ! ----------------------------------------------------------------------------
 
         phi_m = 1.0
@@ -187,7 +186,7 @@ contains
         Pr_t_inv = phi_m / phi_h
 
         ! --- setting output
-        sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+        sfx = sfxDataType(zeta = 0.0, Rib = Rib, &
                 Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
                 Rib_conv_lim = 0.0, &
                 Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
@@ -197,19 +196,18 @@ contains
 
     ! universal functions
     ! --------------------------------------------------------------------------------
-    subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+    subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
         !> @brief universal functions (momentum) & (heat): neutral case
         ! ----------------------------------------------------------------------------
         real, intent(out) :: psi_m, psi_h   !> universal functions
-        real, intent(out) :: zeta           !> = z/L
 
         real, intent(in) :: h0_m, h0_t      !> = z/z0_m, z/z0_h
         real, intent(in) :: B               !> = log(z0_m / z0_h)
         ! ----------------------------------------------------------------------------
 
-        zeta = 0.0
         psi_m = log(h0_m)
         psi_h = log(h0_t) / Pr_t_0_inv
+        !*: this looks redundant z0_t = z0_m in case |B| ~ 0
         if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
 
     end subroutine