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

Last change on this file since 4400 was 4400, checked in by suehring, 15 months ago

Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted

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