source: palm/trunk/SOURCE/basic_constants_and_equations_mod.f90 @ 4811

Last change on this file since 4811 was 4742, checked in by schwenkel, 4 years ago

Implement snow and graupel (bulk microphysics)

  • Property svn:keywords set to Id
File size: 25.7 KB
RevLine 
[3274]1!> @file basic_constants_and_equations_mod.f90
[4509]2!--------------------------------------------------------------------------------------------------!
[3274]3! This file is part of the PALM model system.
4!
[4509]5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
[3274]8!
[4509]9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
[3274]12!
[4509]13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
[3274]15!
[4360]16! Copyright 1997-2020 Leibniz Universitaet Hannover
[4509]17!--------------------------------------------------------------------------------------------------!
[3274]18!
19! Current revisions:
20! -----------------
[4742]21!
22!
[3274]23! Former revisions:
24! -----------------
25! $Id: basic_constants_and_equations_mod.f90 4742 2020-10-14 15:11:02Z maronga $
[4742]26! Implement snow and graupel (bulk microphysics)
27!
28! 4509 2020-04-26 15:57:55Z raasch
[4509]29! file re-formatted to follow the PALM coding standard
30!
31! 4502 2020-04-17 16:14:16Z schwenkel
[4502]32! Implementation of ice microphysics
[4509]33!
[4502]34! 4400 2020-02-10 20:32:41Z suehring
[4509]35! Move routine to transform coordinates from netcdf_interface_mod to
[4400]36! basic_constants_and_equations_mod
[4509]37!
[4400]38! 4360 2020-01-07 11:25:50Z suehring
[4182]39! Corrected "Former revisions" section
[4509]40!
[4182]41! 4088 2019-07-11 13:57:56Z Giersch
[4509]42! Comment of barometric formula improved, function for ideal gas law revised
43!
[4088]44! 4084 2019-07-10 17:09:11Z knoop
[4084]45! Changed precomputed fractions to be variable based
[4509]46!
[4084]47! 4055 2019-06-27 09:47:29Z suehring
[4055]48! Added rgas_univ (universal gas constant) (E.C. Chan)
[4509]49!
50!
[4055]51! 3655 2019-01-07 16:51:22Z knoop
[3634]52! OpenACC port for SPEC
[4182]53! 3361 2018-10-16 20:39:37Z knoop
54! New module (introduced with modularization of bulk cloud physics model)
55!
[3274]56!
[4182]57!
[4509]58!
[3274]59! Description:
60! ------------
[4509]61!> This module contains all basic (physical) constants and functions for the calculation of
62!> diagnostic quantities.
63!-                    -----------------------------------------------------------------------------!
[3274]64 MODULE basic_constants_and_equations_mod
65
66
67    USE kinds
68
69    IMPLICIT NONE
70
[4509]71
[3274]72    REAL(wp), PARAMETER ::  c_p = 1005.0_wp                           !< heat capacity of dry air (J kg-1 K-1)
[4742]73    REAL(wp), PARAMETER ::  c_w = 4185.0_wp                           !< heat capacity of water at 0°C (J kg-1 K-1)
[3449]74    REAL(wp), PARAMETER ::  degc_to_k = 273.15_wp                     !< temperature (in K) of 0 deg C (K)
[3274]75    REAL(wp), PARAMETER ::  g = 9.81_wp                               !< gravitational acceleration (m s-2)
76    REAL(wp), PARAMETER ::  kappa = 0.4_wp                            !< von Karman constant
77    REAL(wp), PARAMETER ::  l_m = 0.33E+06_wp                         !< latent heat of water melting (J kg-1)
78    REAL(wp), PARAMETER ::  l_v = 2.5E+06_wp                          !< latent heat of water vaporization (J kg-1)
79    REAL(wp), PARAMETER ::  l_s = l_m + l_v                           !< latent heat of water sublimation (J kg-1)
80    REAL(wp), PARAMETER ::  molecular_weight_of_nacl = 0.05844_wp     !< mol. m. NaCl (kg mol-1)
81    REAL(wp), PARAMETER ::  molecular_weight_of_c3h4o4 = 0.10406_wp   !< mol. m. malonic acid (kg mol-1)
82    REAL(wp), PARAMETER ::  molecular_weight_of_nh4no3 = 0.08004_wp   !< mol. m. ammonium sulfate (kg mol-1)
83    REAL(wp), PARAMETER ::  molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)
[3449]84    REAL(wp), PARAMETER ::  pi = 3.141592654_wp                       !< PI
[3634]85    !$ACC DECLARE COPYIN(pi)
[4055]86    REAL(wp), PARAMETER ::  rgas_univ = 8.31446261815324_wp           !< universal gas constant (J K-1 mol-1)
[4742]87    REAL(wp), PARAMETER ::  rho_i = 916.7_wp                          !> density of pure ice (kg m-3)
[3274]88    REAL(wp), PARAMETER ::  rho_l = 1.0E3_wp                          !< density of water (kg m-3)
89    REAL(wp), PARAMETER ::  rho_nacl = 2165.0_wp                      !< density of NaCl (kg m-3)
90    REAL(wp), PARAMETER ::  rho_c3h4o4 = 1600.0_wp                    !< density of malonic acid (kg m-3)
91    REAL(wp), PARAMETER ::  rho_nh4no3 = 1720.0_wp                    !< density of ammonium sulfate (kg m-3)
92    REAL(wp), PARAMETER ::  r_d = 287.0_wp                            !< sp. gas const. dry air (J kg-1 K-1)
93    REAL(wp), PARAMETER ::  r_v = 461.51_wp                           !< sp. gas const. water vapor (J kg-1 K-1)
[3449]94    REAL(wp), PARAMETER ::  sigma_sb = 5.67037E-08_wp                 !< Stefan-Boltzmann constant
[3274]95    REAL(wp), PARAMETER ::  solar_constant = 1368.0_wp                !< solar constant at top of atmosphere
96    REAL(wp), PARAMETER ::  vanthoff_nacl = 2.0_wp                    !< van't Hoff factor for NaCl
97    REAL(wp), PARAMETER ::  vanthoff_c3h4o4 = 1.37_wp                 !< van't Hoff factor for malonic acid
98    REAL(wp), PARAMETER ::  vanthoff_nh4no3 = 2.31_wp                 !< van't Hoff factor for ammonium sulfate
99
100    REAL(wp), PARAMETER ::  p_0 = 100000.0_wp                         !< standard pressure reference state
101
[4509]102    REAL(wp), PARAMETER ::  cp_d_rd = c_p / r_d   !< precomputed c_p / r_d
[3274]103    REAL(wp), PARAMETER ::  g_d_cp  = g   / c_p   !< precomputed g / c_p
104    REAL(wp), PARAMETER ::  lv_d_cp = l_v / c_p   !< precomputed l_v / c_p
[4502]105    REAL(wp), PARAMETER ::  ls_d_cp = l_s / c_p   !< precomputed l_s / c_p
[3274]106    REAL(wp), PARAMETER ::  lv_d_rd = l_v / r_d   !< precomputed l_v / r_d
[4084]107    REAL(wp), PARAMETER ::  rd_d_rv = r_d / r_v   !< precomputed r_d / r_v
108    REAL(wp), PARAMETER ::  rd_d_cp = r_d / c_p   !< precomputed r_d / c_p
[3274]109
110    REAL(wp) ::  molecular_weight_of_solute = molecular_weight_of_nacl  !< mol. m. NaCl (kg mol-1)
111    REAL(wp) ::  rho_s = rho_nacl                                       !< density of NaCl (kg m-3)
112    REAL(wp) ::  vanthoff = vanthoff_nacl                               !< van't Hoff factor for NaCl
113
114    SAVE
115
[4509]116    PRIVATE magnus_0d,                                                                             &
117            magnus_1d,                                                                             &
118            magnus_tl_0d,                                                                          &
119            magnus_tl_1d,                                                                          &
120            magnus_0d_ice,                                                                         &
121            magnus_1d_ice,                                                                         &
122            ideal_gas_law_rho_0d,                                                                  &
123            ideal_gas_law_rho_1d,                                                                  &
124            ideal_gas_law_rho_pt_0d,                                                               &
125            ideal_gas_law_rho_pt_1d,                                                               &
126            exner_function_0d,                                                                     &
127            exner_function_1d,                                                                     &
128            exner_function_invers_0d,                                                              &
129            exner_function_invers_1d,                                                              &
130            barometric_formula_0d,                                                                 &
[3274]131            barometric_formula_1d
132
[4400]133
134    INTERFACE convert_utm_to_geographic
135       MODULE PROCEDURE convert_utm_to_geographic
136    END INTERFACE convert_utm_to_geographic
137
[3274]138    INTERFACE magnus
139       MODULE PROCEDURE magnus_0d
140       MODULE PROCEDURE magnus_1d
141    END INTERFACE magnus
142
[4502]143    INTERFACE magnus_tl
144       MODULE PROCEDURE magnus_tl_0d
145       MODULE PROCEDURE magnus_tl_1d
146    END INTERFACE magnus_tl
147
148    INTERFACE magnus_ice
149       MODULE PROCEDURE magnus_0d_ice
150       MODULE PROCEDURE magnus_1d_ice
151    END INTERFACE magnus_ice
152
[3274]153    INTERFACE ideal_gas_law_rho
154       MODULE PROCEDURE ideal_gas_law_rho_0d
155       MODULE PROCEDURE ideal_gas_law_rho_1d
156    END INTERFACE ideal_gas_law_rho
157
158    INTERFACE ideal_gas_law_rho_pt
159       MODULE PROCEDURE ideal_gas_law_rho_pt_0d
160       MODULE PROCEDURE ideal_gas_law_rho_pt_1d
161    END INTERFACE ideal_gas_law_rho_pt
162
163    INTERFACE exner_function
164       MODULE PROCEDURE exner_function_0d
165       MODULE PROCEDURE exner_function_1d
166    END INTERFACE exner_function
167
168    INTERFACE exner_function_invers
169       MODULE PROCEDURE exner_function_invers_0d
170       MODULE PROCEDURE exner_function_invers_1d
171    END INTERFACE exner_function_invers
172
173    INTERFACE barometric_formula
174       MODULE PROCEDURE barometric_formula_0d
175       MODULE PROCEDURE barometric_formula_1d
176    END INTERFACE barometric_formula
[4400]177!
178!-- Public routines
179    PUBLIC convert_utm_to_geographic
[3274]180
181 CONTAINS
182
[4400]183
[4509]184!--------------------------------------------------------------------------------------------------!
[3274]185! Description:
186! ------------
[4509]187!> Convert UTM coordinates into geographic latitude and longitude. Conversion is based on the work
188!> of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28 and Karney (2013) DOI: 10.1007/s00190-012-0578-z
[4400]189!> Based on a JavaScript of the geodesy function library written by chrisveness
190!> https://github.com/chrisveness/geodesy
[4509]191!--------------------------------------------------------------------------------------------------!
192 SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
[4400]193
[4509]194    INTEGER(iwp) ::  j   !< loop index
[4400]195
[4509]196    REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
197    REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
198    REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
199    REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
[4400]200
[4509]201    REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
202    REAL(wp) ::  cos_eta_s   !< cos(eta_s)
203    REAL(wp) ::  delta_i     !<
204    REAL(wp) ::  delta_tau_i !<
205    REAL(wp) ::  e           !< eccentricity
206    REAL(wp) ::  eta         !<
207    REAL(wp) ::  eta_s       !<
208    REAL(wp) ::  n           !< 3rd flattening
209    REAL(wp) ::  n2          !< n^2
210    REAL(wp) ::  n3          !< n^3
211    REAL(wp) ::  n4          !< n^4
212    REAL(wp) ::  n5          !< n^5
213    REAL(wp) ::  n6          !< n^6
214    REAL(wp) ::  nu          !<
215    REAL(wp) ::  nu_s        !<
216    REAL(wp) ::  sin_eta_s   !< sin(eta_s)
217    REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
218    REAL(wp) ::  tau_i       !<
219    REAL(wp) ::  tau_i_s     !<
220    REAL(wp) ::  tau_s       !<
221    REAL(wp) ::  x           !< adjusted easting
222    REAL(wp) ::  y           !< adjusted northing
[4400]223
[4509]224    REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
[4400]225
[4509]226    REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
227                                               !< (/semi_major_axis,
228                                               !<   inverse_flattening,
229                                               !<   longitude_of_prime_meridian,
230                                               !<   longitude_of_central_meridian,
231                                               !<   scale_factor_at_central_meridian,
232                                               !<   latitude_of_projection_origin,
233                                               !<   false_easting,
234                                               !<   false_northing /)
[4400]235
[4509]236    x = eutm - crs(7)  ! remove false easting
237    y = nutm - crs(8)  ! remove false northing
[4400]238!
[4509]239!-- From Karney 2011 Eq 15-22, 36:
240    e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
241    n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
242    n2 = n * n
243    n3 = n * n2
244    n4 = n * n3
245    n5 = n * n4
246    n6 = n * n5
[4400]247
[4509]248    a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2 + 0.015625_wp * n4 + 3.90625E-3_wp * n6 )
[4400]249
[4509]250    nu  = x / ( crs(5) * a )
251    eta = y / ( crs(5) * a )
[4400]252
[4509]253!-- According to KrÃŒger (1912), eq. 26*
254    beta(1) =          0.5_wp                  * n                                                 &
255              -        2.0_wp /         3.0_wp * n2                                                &
256              +       37.0_wp /        96.0_wp * n3                                                &
257              -        1.0_wp /       360.0_wp * n4                                                &
258              -       81.0_wp /       512.0_wp * n5                                                &
259              +    96199.0_wp /    604800.0_wp * n6
[4400]260
[4509]261    beta(2) =          1.0_wp /        48.0_wp * n2                                                &
262              +        1.0_wp /        15.0_wp * n3                                                &
263              -      437.0_wp /      1440.0_wp * n4                                                &
264              +       46.0_wp /       105.0_wp * n5                                                &
265              -  1118711.0_wp /   3870720.0_wp * n6
[4400]266
[4509]267    beta(3) =         17.0_wp /       480.0_wp * n3                                                &
268              -       37.0_wp /       840.0_wp * n4                                                &
269              -      209.0_wp /      4480.0_wp * n5                                                &
270              +     5569.0_wp /     90720.0_wp * n6
[4400]271
[4509]272    beta(4) =       4397.0_wp /    161280.0_wp * n4                                                &
273              -       11.0_wp /       504.0_wp * n5                                                &
274              -   830251.0_wp /   7257600.0_wp * n6
[4400]275
[4509]276    beta(5) =       4583.0_wp /    161280.0_wp * n5                                                &
277              -   108847.0_wp /   3991680.0_wp * n6
[4400]278
[4509]279    beta(6) =   20648693.0_wp / 638668800.0_wp * n6
[4400]280
[4509]281    eta_s = eta
282    nu_s  = nu
283    DO  j = 1, 6
284      eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
285      nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
286    ENDDO
[4400]287
[4509]288    sinh_nu_s = SINH( nu_s )
289    sin_eta_s = SIN( eta_s )
290    cos_eta_s = COS( eta_s )
[4400]291
[4509]292    tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
[4400]293
[4509]294    tau_i = tau_s
295    delta_tau_i = 1.0_wp
[4400]296
[4509]297    DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
[4400]298
[4509]299      delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
[4400]300
[4509]301      tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 ) - delta_i * SQRT( 1.0_wp + tau_i**2 )
[4400]302
[4509]303      delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )                              &
304                    * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )                                    &
305                    / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
[4400]306
[4509]307      tau_i = tau_i + delta_tau_i
[4400]308
[4509]309    ENDDO
[4400]310
[4509]311    lat = ATAN( tau_i ) / pi * 180.0_wp
312    lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
[4400]313
[4509]314 END SUBROUTINE convert_utm_to_geographic
[4400]315
[4509]316!--------------------------------------------------------------------------------------------------!
[4400]317! Description:
318! ------------
[3274]319!> This function computes the magnus formula (Press et al., 1992).
[4509]320!> The magnus formula is needed to calculate the saturation vapor pressure.
321!--------------------------------------------------------------------------------------------------!
322 FUNCTION magnus_0d( t )
[3274]323
[4509]324    IMPLICIT NONE
[3274]325
[4509]326    REAL(wp), INTENT(IN) ::  t  !< temperature (K)
[3274]327
[4509]328    REAL(wp) ::  magnus_0d
329
[3274]330!
[4509]331!-- Saturation vapor pressure for a specific temperature:
332    magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) / ( t - 29.65_wp  ) )
[3274]333
[4509]334 END FUNCTION magnus_0d
[3274]335
[4509]336!--------------------------------------------------------------------------------------------------!
[3274]337! Description:
338! ------------
339!> This function computes the magnus formula (Press et al., 1992).
[4509]340!> The magnus formula is needed to calculate the saturation vapor pressure.
341!--------------------------------------------------------------------------------------------------!
342 FUNCTION magnus_1d( t )
[3274]343
[4509]344    IMPLICIT NONE
[3274]345
[4509]346    REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
[3274]347
[4509]348    REAL(wp), DIMENSION(size(t)) ::  magnus_1d
349
[3274]350!
[4509]351!-- Saturation vapor pressure for a specific temperature:
352    magnus_1d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) / ( t - 29.65_wp  ) )
[3274]353
[4509]354 END FUNCTION magnus_1d
[3274]355
[4509]356!--------------------------------------------------------------------------------------------------!
[3274]357! Description:
358! ------------
[4509]359!> This function computes the magnus formula (Press et al., 1992) using the (ice-) liquid water
360!> potential temperature.
361!> The magnus formula is needed to calculate the saturation vapor pressure over a plane liquid water
362!> surface.
363!--------------------------------------------------------------------------------------------------!
364 FUNCTION magnus_tl_0d( t_l )
[4502]365
[4509]366    IMPLICIT NONE
[4502]367
[4509]368    REAL(wp), INTENT(IN) ::  t_l  !< liquid water temperature (K)
[4502]369
[4509]370    REAL(wp) ::  magnus_tl_0d
371
[4502]372!
[4509]373!-- Saturation vapor pressure for a specific temperature:
374    magnus_tl_0d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp  ) )
[4502]375
[4509]376 END FUNCTION magnus_tl_0d
[4502]377
[4509]378!--------------------------------------------------------------------------------------------------!
[4502]379! Description:
380! ------------
[4509]381!> This function computes the magnus formula (Press et al., 1992) using the (ice-) liquid water
382!> potential temperature.
383!> The magnus formula is needed to calculate the saturation vapor pressure over a plane liquid water
384!> surface.
385!--------------------------------------------------------------------------------------------------!
386 FUNCTION magnus_tl_1d( t_l )
[4502]387
[4509]388    IMPLICIT NONE
[4502]389
[4509]390    REAL(wp), INTENT(IN), DIMENSION(:) ::  t_l  !< liquid water temperature (K)
[4502]391
[4509]392    REAL(wp), DIMENSION(size(t_l)) ::  magnus_tl_1d
[4502]393!
[4509]394!-- Saturation vapor pressure for a specific temperature:
395    magnus_tl_1d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp  ) )
[4502]396
[4509]397 END FUNCTION magnus_tl_1d
[4502]398
[4509]399!--------------------------------------------------------------------------------------------------!
[4502]400! Description:
401! ------------
402!> This function computes the magnus formula (Press et al., 1992).
[4509]403!> The magnus formula is needed to calculate the saturation vapor pressure over a plane ice surface.
404!--------------------------------------------------------------------------------------------------!
405 FUNCTION magnus_0d_ice( t )
[4502]406
[4509]407    IMPLICIT NONE
[4502]408
[4509]409    REAL(wp), INTENT(IN) ::  t  !< temperature (K)
[4502]410
[4509]411    REAL(wp) ::  magnus_0d_ice
412
[4502]413!
[4742]414!-- Saturation vapor pressure for a specific temperature:
415    !magnus_0d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) / ( t - 0.53_wp  ) )
416    magnus_0d_ice =  610.78_wp * EXP( 21.875_wp * ( t - degc_to_k ) / ( t - 7.66_wp  ) )
[4502]417
[4742]418
[4509]419 END FUNCTION magnus_0d_ice
[4502]420
[4509]421!--------------------------------------------------------------------------------------------------!
[4502]422! Description:
423! ------------
424!> This function computes the magnus formula (Press et al., 1992).
[4509]425!> The magnus formula is needed to calculate the saturation vapor pressure over a plane ice surface.
426!--------------------------------------------------------------------------------------------------!
427 FUNCTION magnus_1d_ice( t )
[4502]428
[4509]429    IMPLICIT NONE
[4502]430
[4509]431    REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
[4502]432
[4509]433    REAL(wp), DIMENSION(size(t)) ::  magnus_1d_ice
434
[4502]435!
[4509]436!-- Saturation vapor pressure for a specific temperature:
[4742]437    !magnus_1d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) / ( t - 0.53_wp  ) )
438    magnus_1d_ice =  610.78_wp * EXP( 21.875_wp * ( t - degc_to_k ) / ( t - 7.66_wp  ) )
[4502]439
[4742]440
[4509]441 END FUNCTION magnus_1d_ice
[4502]442
[4509]443!--------------------------------------------------------------------------------------------------!
[4502]444! Description:
445! ------------
[3274]446!> Compute the ideal gas law for scalar arguments.
[4509]447!--------------------------------------------------------------------------------------------------!
448 FUNCTION ideal_gas_law_rho_0d( p, t )
[3274]449
[4509]450    IMPLICIT NONE
[3274]451
[4509]452    REAL(wp), INTENT(IN) ::  p  !< pressure (Pa)
453    REAL(wp), INTENT(IN) ::  t  !< temperature (K)
[3274]454
[4509]455    REAL(wp) ::  ideal_gas_law_rho_0d
456
[3274]457!
[4509]458!-- Compute density according to ideal gas law:
459    ideal_gas_law_rho_0d = p / (r_d * t)
[3274]460
[4509]461 END FUNCTION ideal_gas_law_rho_0d
[3274]462
[4509]463!--------------------------------------------------------------------------------------------------!
[3274]464! Description:
465! ------------
466!> Compute the ideal gas law for 1-D array arguments.
[4509]467!--------------------------------------------------------------------------------------------------!
468 FUNCTION ideal_gas_law_rho_1d( p, t )
[3274]469
[4509]470    IMPLICIT NONE
[3274]471
[4509]472    REAL(wp), INTENT(IN), DIMENSION(:) ::  p  !< pressure (Pa)
473    REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
[3274]474
[4509]475    REAL(wp), DIMENSION(size(p)) ::  ideal_gas_law_rho_1d
476
[3274]477!
[4509]478!-- Compute density according to ideal gas law:
479    ideal_gas_law_rho_1d = p / (r_d * t)
[3274]480
[4509]481 END FUNCTION ideal_gas_law_rho_1d
[3274]482
[4509]483!--------------------------------------------------------------------------------------------------!
[3274]484! Description:
485! ------------
486!> Compute the ideal gas law for scalar arguments.
[4509]487!--------------------------------------------------------------------------------------------------!
488 FUNCTION ideal_gas_law_rho_pt_0d( p, t )
[3274]489
[4509]490    IMPLICIT NONE
[3274]491
[4509]492    REAL(wp), INTENT(IN) ::  p  !< pressure (Pa)
493    REAL(wp), INTENT(IN) ::  t  !< temperature (K)
[3274]494
[4509]495    REAL(wp) ::  ideal_gas_law_rho_pt_0d
496
[3274]497!
[4509]498!-- Compute density according to ideal gas law:
499    ideal_gas_law_rho_pt_0d = p / (r_d * exner_function(p) * t)
[3274]500
[4509]501 END FUNCTION ideal_gas_law_rho_pt_0d
[3274]502
[4509]503!--------------------------------------------------------------------------------------------------!
[3274]504! Description:
505! ------------
506!> Compute the ideal gas law for 1-D array arguments.
[4509]507!--------------------------------------------------------------------------------------------------!
508 FUNCTION ideal_gas_law_rho_pt_1d( p, t )
[3274]509
[4509]510    IMPLICIT NONE
[3274]511
[4509]512    REAL(wp), INTENT(IN), DIMENSION(:) ::  p  !< pressure (Pa)
513    REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
[3274]514
[4509]515    REAL(wp), DIMENSION(size(p)) ::  ideal_gas_law_rho_pt_1d
516
[3274]517!
[4509]518!-- Compute density according to ideal gas law:
519    ideal_gas_law_rho_pt_1d = p / (r_d * exner_function(p) * t)
[3274]520
[4509]521 END FUNCTION ideal_gas_law_rho_pt_1d
[3274]522
[4509]523!--------------------------------------------------------------------------------------------------!
[3274]524! Description:
525! ------------
526!> Compute the exner function for scalar arguments.
[4509]527!--------------------------------------------------------------------------------------------------!
528 FUNCTION exner_function_0d( p )
[3274]529
[4509]530    IMPLICIT NONE
[3274]531
[4509]532    REAL(wp), INTENT(IN) ::  p    !< pressure (Pa)
[3274]533
[4509]534    REAL(wp) ::  exner_function_0d
535
[3274]536!
[4509]537!-- Compute exner function:
538    exner_function_0d = ( p / p_0 )**( rd_d_cp )
[3274]539
[4509]540 END FUNCTION exner_function_0d
[3274]541
[4509]542!--------------------------------------------------------------------------------------------------!
[3274]543! Description:
544! ------------
545!> Compute the exner function for 1-D array arguments.
[4509]546!--------------------------------------------------------------------------------------------------!
547 FUNCTION exner_function_1d( p )
[3274]548
[4509]549    IMPLICIT NONE
[3274]550
[4509]551    REAL(wp), INTENT(IN), DIMENSION(:) ::  p  !< pressure (Pa)
[3274]552
[4509]553    REAL(wp), DIMENSION(size(p)) ::  exner_function_1d
554
[3274]555!
[4509]556!-- Compute exner function:
557    exner_function_1d = ( p / p_0 )**( rd_d_cp )
[3274]558
[4509]559 END FUNCTION exner_function_1d
[3274]560
[4509]561!--------------------------------------------------------------------------------------------------!
[3274]562! Description:
563! ------------
564!> Compute the exner function for scalar arguments.
[4509]565!--------------------------------------------------------------------------------------------------!
566 FUNCTION exner_function_invers_0d( p )
[3274]567
[4509]568    IMPLICIT NONE
[3274]569
[4509]570    REAL(wp), INTENT(IN) ::  p    !< pressure (Pa)
[3274]571
[4509]572    REAL(wp) ::  exner_function_invers_0d
573
[3274]574!
[4509]575!-- Compute exner function:
576    exner_function_invers_0d = ( p_0 / p )**( rd_d_cp )
[3274]577
[4509]578 END FUNCTION exner_function_invers_0d
[3274]579
[4509]580!--------------------------------------------------------------------------------------------------!
[3274]581! Description:
582! ------------
583!> Compute the exner function for 1-D array arguments.
[4509]584!--------------------------------------------------------------------------------------------------!
585 FUNCTION exner_function_invers_1d( p )
[3274]586
[4509]587    IMPLICIT NONE
[3274]588
[4509]589    REAL(wp), INTENT(IN), DIMENSION(:) ::  p  !< pressure (Pa)
[3274]590
[4509]591    REAL(wp), DIMENSION(size(p)) ::  exner_function_invers_1d
592
[3274]593!
[4509]594!-- Compute exner function:
595    exner_function_invers_1d = ( p_0 / p )**( rd_d_cp )
[3274]596
[4509]597 END FUNCTION exner_function_invers_1d
[3274]598
[4509]599!--------------------------------------------------------------------------------------------------!
[3274]600! Description:
601! ------------
[4509]602!> Compute the barometric formula for scalar arguments. The calculation is based on the assumption
603!> of a polytropic atmosphere and neutral stratification, where the temperature lapse rate is g/cp.
604!--------------------------------------------------------------------------------------------------!
605 FUNCTION barometric_formula_0d( z, t_0, p_0)
[3274]606
[4509]607    IMPLICIT NONE
[3274]608
[4509]609    REAL(wp), INTENT(IN) ::  z    !< height (m)
610    REAL(wp), INTENT(IN) ::  t_0  !< temperature reference state (K)
611    REAL(wp), INTENT(IN) ::  p_0  !< surface pressure (Pa)
[3274]612
[4509]613    REAL(wp) ::  barometric_formula_0d
614
[3274]615!
[4509]616!-- Compute barometric formula:
617    barometric_formula_0d =  p_0 * ( (t_0 - g_d_cp * z) / t_0 )**( cp_d_rd )
[3274]618
[4509]619 END FUNCTION barometric_formula_0d
[3274]620
[4509]621!--------------------------------------------------------------------------------------------------!
[3274]622! Description:
623! ------------
[4509]624!> Compute the barometric formula for 1-D array arguments. The calculation is based on the
625!> assumption of a polytropic atmosphere and neutral stratification, where the temperature lapse
626!> rate is g/cp.
627!--------------------------------------------------------------------------------------------------!
628 FUNCTION barometric_formula_1d( z, t_0, p_0)
[3274]629
[4509]630    IMPLICIT NONE
[3274]631
[4509]632    REAL(wp), INTENT(IN), DIMENSION(:) ::  z  !< height (m)
633    REAL(wp), INTENT(IN) ::  t_0              !< temperature reference state (K)
634    REAL(wp), INTENT(IN) ::  p_0              !< surface pressure (Pa)
[3274]635
[4509]636    REAL(wp), DIMENSION(size(z)) ::  barometric_formula_1d
637
[3274]638!
[4509]639!-- Compute barometric formula:
640    barometric_formula_1d =  p_0 * ( (t_0 - g_d_cp * z) / t_0 )**( cp_d_rd )
[3274]641
[4509]642 END FUNCTION barometric_formula_1d
[3274]643
644 END MODULE basic_constants_and_equations_mod
Note: See TracBrowser for help on using the repository browser.