Notes on the implementation in PALM
The following code is from PALM
❯ svn info
Path: .
Working Copy Root Path: /home/avmo/src/snek5000/misc/palm
URL: https://palm.muk.uni-hannover.de/svn/palm/trunk
Relative URL: ^/palm/trunk
Repository Root: https://palm.muk.uni-hannover.de/svn
Repository UUID: 5fb83930-4629-0410-b46d-57b98ee623a8
Revision: 4901
Node Kind: directory
Schedule: normal
Last Changed Author: banzhafs
Last Changed Rev: 4901
Last Changed Date: 2021-03-04 22:24:08 +0100 (Thu, 04 Mar 2021)
This excerpt from subroutine calc_ob is responsible for computing the Obukhov Length.
!> @file surface_layer_fluxes_mod.f90
!--------------------------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
! Public License as published by the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
! Public License for more details.
!
! You should have received a copy of the GNU General Public License along with PALM. If not, see
! <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2021 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------------------------!
!-- Calculate the Obukhov length using Newton iteration
!$OMP PARALLEL DO PRIVATE(i, j, z_mo, ol_old, iter, ol_m, ol_l, ol_u, f, f_d_ol)
!$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
!$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
!$ACC PRESENT(surf)
DO m = 1, surf%ns
i = surf%i(m)
j = surf%j(m)
z_mo = surf%z_mo(m)
!
!-- Store current value in case the Newton iteration fails
ol_old = surf%ol(m)
!
!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign
IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN
IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp
IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
ENDIF
!
!-- Iteration to find Obukhov length
iter = 0
DO
iter = iter + 1
!
!-- In case of divergence, use the value of the previous time step
IF ( iter > 1000 ) THEN
surf%ol(m) = ol_old
EXIT
ENDIF
ol_m = surf%ol(m)
ol_l = ol_m - 0.001_wp * ol_m
ol_u = ol_m + 0.001_wp * ol_m
IF ( ibc_pt_b /= 1 ) THEN
!
!-- Calculate f = Ri - [...]/[...]^2 = 0
f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_z0h(m) &
- psi_h( z_mo / ol_m ) &
+ psi_h( surf%z0h(m) / ol_m ) ) / &
( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) &
+ psi_m( surf%z0(m) / ol_m ) )**2
!
!-- Calculate df/dL
f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_z0h(m) &
- psi_h( z_mo / ol_u ) &
+ psi_h( surf%z0h(m) / ol_u ) ) / &
( surf%ln_z_z0(m) - psi_m( z_mo / ol_u ) &
+ psi_m( surf%z0(m) / ol_u ) )**2 &
+ ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l ) &
+ psi_h( surf%z0h(m) / ol_l ) ) /&
( surf%ln_z_z0(m) - psi_m( z_mo / ol_l ) &
+ psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )
ELSE
!
!-- Calculate f = Ri - 1 /[...]^3 = 0
f = surf%rib(m) - ( z_mo / ol_m ) / &
( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) )**3
!
!-- Calculate df/dL
f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_z0(m) &
- psi_m( z_mo / ol_u ) &
+ psi_m( surf%z0(m) / ol_u ) &
)**3 &
+ ( z_mo / ol_l ) / ( surf%ln_z_z0(m) &
- psi_m( z_mo / ol_l ) &
+ psi_m( surf%z0(m) / ol_l ) &
)**3 &
) / ( ol_u - ol_l )
ENDIF
!
!-- Calculate new L
surf%ol(m) = ol_m - f / f_d_ol
!
!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and
!-- ensure convergence.
IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp
!
!-- If unrealistic value occurs, set L to the maximum value that is allowed
IF ( ABS( surf%ol(m) ) > ol_max ) THEN
surf%ol(m) = ol_max
EXIT
ENDIF
!
!-- Assure that Obukhov length does not become zero. If the limit is reached, exit the loop.
IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN
surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
EXIT
ENDIF
!
!-- Check for convergence
IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT
ENDDO
ENDDO
What is interesting here is that f_d_ol is not an exact derivative but, uses approximates it using Central difference. Does not seem like Secant Method.