{ "cells": [ { "cell_type": "markdown", "id": "f6096238-fdfe-4f91-8d25-89279494b577", "metadata": {}, "source": [ "# Notes on the implementation in PALM\n", "\n", "The following code is from PALM\n", "```\n", "❯ svn info\n", "Path: .\n", "Working Copy Root Path: /home/avmo/src/snek5000/misc/palm\n", "URL: https://palm.muk.uni-hannover.de/svn/palm/trunk\n", "Relative URL: ^/palm/trunk\n", "Repository Root: https://palm.muk.uni-hannover.de/svn\n", "Repository UUID: 5fb83930-4629-0410-b46d-57b98ee623a8\n", "Revision: 4901\n", "Node Kind: directory\n", "Schedule: normal\n", "Last Changed Author: banzhafs\n", "Last Changed Rev: 4901\n", "Last Changed Date: 2021-03-04 22:24:08 +0100 (Thu, 04 Mar 2021)\n", "\n", "```\n", "\n", "This excerpt from subroutine `calc_ob` is responsible for computing the Obukhov Length.\n", "\n", "```fortran\n", "!> @file surface_layer_fluxes_mod.f90\n", "!--------------------------------------------------------------------------------------------------!\n", "! This file is part of the PALM model system.\n", "!\n", "! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General\n", "! Public License as published by the Free Software Foundation, either version 3 of the License, or\n", "! (at your option) any later version.\n", "!\n", "! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the\n", "! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General\n", "! Public License for more details.\n", "!\n", "! You should have received a copy of the GNU General Public License along with PALM. If not, see\n", "! .\n", "!\n", "! Copyright 1997-2021 Leibniz Universitaet Hannover\n", "!--------------------------------------------------------------------------------------------------!\n", "\n", "\n", "!-- Calculate the Obukhov length using Newton iteration\n", " !$OMP PARALLEL DO PRIVATE(i, j, z_mo, ol_old, iter, ol_m, ol_l, ol_u, f, f_d_ol)\n", " !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &\n", " !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &\n", " !$ACC PRESENT(surf)\n", " DO m = 1, surf%ns\n", " i = surf%i(m)\n", " j = surf%j(m)\n", "\n", " z_mo = surf%z_mo(m)\n", "\n", "!\n", "!-- Store current value in case the Newton iteration fails\n", " ol_old = surf%ol(m)\n", "\n", "!\n", "!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign\n", " IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN\n", " IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp\n", " IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp\n", " ENDIF\n", "!\n", "!-- Iteration to find Obukhov length\n", " iter = 0\n", " DO\n", " iter = iter + 1\n", "!\n", "!-- In case of divergence, use the value of the previous time step\n", " IF ( iter > 1000 ) THEN\n", " surf%ol(m) = ol_old\n", " EXIT\n", " ENDIF\n", "\n", " ol_m = surf%ol(m)\n", " ol_l = ol_m - 0.001_wp * ol_m\n", " ol_u = ol_m + 0.001_wp * ol_m\n", "\n", "\n", " IF ( ibc_pt_b /= 1 ) THEN\n", "!\n", "!-- Calculate f = Ri - [...]/[...]^2 = 0\n", " f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_z0h(m) &\n", " - psi_h( z_mo / ol_m ) &\n", " + psi_h( surf%z0h(m) / ol_m ) ) / &\n", " ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) &\n", " + psi_m( surf%z0(m) / ol_m ) )**2\n", "\n", "!\n", "!-- Calculate df/dL\n", " f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_z0h(m) &\n", " - psi_h( z_mo / ol_u ) &\n", " + psi_h( surf%z0h(m) / ol_u ) ) / &\n", " ( surf%ln_z_z0(m) - psi_m( z_mo / ol_u ) &\n", " + psi_m( surf%z0(m) / ol_u ) )**2 &\n", " + ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l ) &\n", " + psi_h( surf%z0h(m) / ol_l ) ) /&\n", " ( surf%ln_z_z0(m) - psi_m( z_mo / ol_l ) &\n", " + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )\n", " ELSE\n", "!\n", "!-- Calculate f = Ri - 1 /[...]^3 = 0\n", " f = surf%rib(m) - ( z_mo / ol_m ) / &\n", " ( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) )**3\n", "\n", "!\n", "!-- Calculate df/dL\n", " f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_z0(m) &\n", " - psi_m( z_mo / ol_u ) &\n", " + psi_m( surf%z0(m) / ol_u ) &\n", " )**3 &\n", " + ( z_mo / ol_l ) / ( surf%ln_z_z0(m) &\n", " - psi_m( z_mo / ol_l ) &\n", " + psi_m( surf%z0(m) / ol_l ) &\n", " )**3 &\n", " ) / ( ol_u - ol_l )\n", " ENDIF\n", "!\n", "!-- Calculate new L\n", " surf%ol(m) = ol_m - f / f_d_ol\n", "\n", "!\n", "!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and\n", "!-- ensure convergence.\n", " IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp\n", "!\n", "!-- If unrealistic value occurs, set L to the maximum value that is allowed\n", " IF ( ABS( surf%ol(m) ) > ol_max ) THEN\n", " surf%ol(m) = ol_max\n", " EXIT\n", " ENDIF\n", "!\n", "!-- Assure that Obukhov length does not become zero. If the limit is reached, exit the loop.\n", " IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN\n", " surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )\n", " EXIT\n", " ENDIF\n", "!\n", "!-- Check for convergence\n", " IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT\n", " ENDDO\n", " ENDDO\n", "```" ] }, { "cell_type": "markdown", "id": "a2f889e8-2f11-47fb-8717-b8fa2427a4b6", "metadata": {}, "source": [ "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](https://jaketae.github.io/study/numerical-methods/). " ] }, { "cell_type": "code", "execution_count": null, "id": "294fd5bd-82ec-4fb7-9132-728057f7602a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "py-snek", "language": "python", "name": "py-snek" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.9" } }, "nbformat": 4, "nbformat_minor": 5 }