Changeset 4400 for palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

Timestamp:

Feb 10, 2020 8:32:41 PM (4 years ago)

Author:

suehring

Message:

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

File:

• palm/trunk/SOURCE/basic_constants_and_equations_mod.f90 (modified) (4 diffs)

palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Move routine to transform coordinates from netcdf_interface_mod to 
+! basic_constants_and_equations_mod
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -113 +117 @@
             barometric_formula_1d
 
+
+    INTERFACE convert_utm_to_geographic
+       MODULE PROCEDURE convert_utm_to_geographic
+    END INTERFACE convert_utm_to_geographic
+
     INTERFACE magnus
        MODULE PROCEDURE magnus_0d
@@ -142 +151 @@
        MODULE PROCEDURE barometric_formula_1d
     END INTERFACE barometric_formula
+!
+!-- Public routines
+    PUBLIC convert_utm_to_geographic
 
  CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Convert UTM coordinates into geographic latitude and longitude. Conversion
+!> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
+!> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
+!> Based on a JavaScript of the geodesy function library written by chrisveness
+!> https://github.com/chrisveness/geodesy
+!------------------------------------------------------------------------------!
+    SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
+
+       INTEGER(iwp) ::  j   !< loop index
+
+       REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
+       REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
+       REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
+       REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
+
+       REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
+       REAL(wp) ::  cos_eta_s   !< cos(eta_s)
+       REAL(wp) ::  delta_i     !<
+       REAL(wp) ::  delta_tau_i !<
+       REAL(wp) ::  e           !< eccentricity
+       REAL(wp) ::  eta         !<
+       REAL(wp) ::  eta_s       !<
+       REAL(wp) ::  n           !< 3rd flattening
+       REAL(wp) ::  n2          !< n^2
+       REAL(wp) ::  n3          !< n^3
+       REAL(wp) ::  n4          !< n^4
+       REAL(wp) ::  n5          !< n^5
+       REAL(wp) ::  n6          !< n^6
+       REAL(wp) ::  nu          !<
+       REAL(wp) ::  nu_s        !<
+       REAL(wp) ::  sin_eta_s   !< sin(eta_s)
+       REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
+       REAL(wp) ::  tau_i       !<
+       REAL(wp) ::  tau_i_s     !<
+       REAL(wp) ::  tau_s       !<
+       REAL(wp) ::  x           !< adjusted easting
+       REAL(wp) ::  y           !< adjusted northing
+
+       REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
+
+       REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
+                                                  !< (/semi_major_axis,
+                                                  !<   inverse_flattening,
+                                                  !<   longitude_of_prime_meridian,
+                                                  !<   longitude_of_central_meridian,
+                                                  !<   scale_factor_at_central_meridian,
+                                                  !<   latitude_of_projection_origin,
+                                                  !<   false_easting,
+                                                  !<   false_northing /)
+
+       x = eutm - crs(7)  ! remove false easting
+       y = nutm - crs(8)  ! remove false northing
+!
+!--    from Karney 2011 Eq 15-22, 36:
+       e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
+       n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
+       n2 = n * n
+       n3 = n * n2
+       n4 = n * n3
+       n5 = n * n4
+       n6 = n * n5
+
+       a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
+                                              + 0.015625_wp * n4   &
+                                              + 3.90625E-3_wp * n6 )
+
+       nu  = x / ( crs(5) * a )
+       eta = y / ( crs(5) * a )
+
+!--    According to KrÃŒger (1912), eq. 26*
+       beta(1) =        0.5_wp                  * n  &
+               -        2.0_wp /         3.0_wp * n2 &
+               +       37.0_wp /        96.0_wp * n3 &
+               -        1.0_wp /       360.0_wp * n4 &
+               -       81.0_wp /       512.0_wp * n5 &
+               +    96199.0_wp /    604800.0_wp * n6
+
+       beta(2) =        1.0_wp /        48.0_wp * n2 &
+               +        1.0_wp /        15.0_wp * n3 &
+               -      437.0_wp /      1440.0_wp * n4 &
+               +       46.0_wp /       105.0_wp * n5 &
+               -  1118711.0_wp /   3870720.0_wp * n6
+
+       beta(3) =       17.0_wp /       480.0_wp * n3 &
+               -       37.0_wp /       840.0_wp * n4 &
+               -      209.0_wp /      4480.0_wp * n5 &
+               +     5569.0_wp /     90720.0_wp * n6
+
+       beta(4) =     4397.0_wp /    161280.0_wp * n4 &
+               -       11.0_wp /       504.0_wp * n5 &
+               -   830251.0_wp /   7257600.0_wp * n6
+
+       beta(5) =     4583.0_wp /    161280.0_wp * n5 &
+               -   108847.0_wp /   3991680.0_wp * n6
+
+       beta(6) = 20648693.0_wp / 638668800.0_wp * n6
+
+       eta_s = eta
+       nu_s  = nu
+       DO  j = 1, 6
+         eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
+         nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
+       ENDDO
+
+       sinh_nu_s = SINH( nu_s )
+       sin_eta_s = SIN( eta_s )
+       cos_eta_s = COS( eta_s )
+
+       tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
+
+       tau_i = tau_s
+       delta_tau_i = 1.0_wp
+
+       DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
+
+         delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
+
+         tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
+                  - delta_i * SQRT( 1.0_wp + tau_i**2 )
+
+         delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
+                      * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
+                      / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
+
+         tau_i = tau_i + delta_tau_i
+
+       ENDDO
+
+       lat = ATAN( tau_i ) / pi * 180.0_wp
+       lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
+
+    END SUBROUTINE convert_utm_to_geographic
 
 !------------------------------------------------------------------------------!
@@ -161 +310 @@
 !--    Saturation vapor pressure for a specific temperature:
        magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) /             &
-                                                   ( t - 29.65_wp  ) )
+                                               ( t - 29.65_wp  ) )
 
     END FUNCTION magnus_0d