Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r4360 r4400  
    2525! -----------------
    2626! $Id$
     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
    2731! Adjusted output of multi-agent system for biometeorology
    2832!
     
    135139
    136140    USE netcdf_data_input_mod,                                                 &
    137         ONLY: coord_ref_sys, init_model
     141        ONLY: coord_ref_sys,                                                   &
     142              crs_list,                                                        &
     143              init_model
    138144
    139145    PRIVATE
     
    417423
    418424    USE basic_constants_and_equations_mod,                                     &
    419         ONLY:  pi
     425        ONLY:  convert_utm_to_geographic,                                      &
     426               pi
    420427
    421428    USE control_parameters,                                                    &
     
    564571
    565572    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
    566     REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
    567573
    568574    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
     
    662668
    663669    ENDIF
    664 !
    665 !-- Convert coord_ref_sys into vector (used for lat/lon calculation)
    666     crs_list = (/ coord_ref_sys%semi_major_axis,                  &
    667                   coord_ref_sys%inverse_flattening,               &
    668                   coord_ref_sys%longitude_of_prime_meridian,      &
    669                   coord_ref_sys%longitude_of_central_meridian,    &
    670                   coord_ref_sys%scale_factor_at_central_meridian, &
    671                   coord_ref_sys%latitude_of_projection_origin,    &
    672                   coord_ref_sys%false_easting,                    &
    673                   coord_ref_sys%false_northing /)
    674670
    675671!
     
    71087104
    71097105
    7110 !------------------------------------------------------------------------------!
    7111 ! Description:
    7112 ! ------------
    7113 !> Convert UTM coordinates into geographic latitude and longitude. Conversion
    7114 !> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
    7115 !> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
    7116 !> Based on a JavaScript of the geodesy function library written by chrisveness
    7117 !> https://github.com/chrisveness/geodesy
    7118 !------------------------------------------------------------------------------!
    7119  SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
    7120 
    7121     USE basic_constants_and_equations_mod,                                     &
    7122         ONLY:  pi
    7123 
    7124     IMPLICIT NONE
    7125 
    7126     INTEGER(iwp) ::  j   !< loop index
    7127 
    7128     REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
    7129     REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
    7130     REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
    7131     REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
    7132 
    7133     REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
    7134     REAL(wp) ::  cos_eta_s   !< cos(eta_s)
    7135     REAL(wp) ::  delta_i     !<
    7136     REAL(wp) ::  delta_tau_i !<
    7137     REAL(wp) ::  e           !< eccentricity
    7138     REAL(wp) ::  eta         !<
    7139     REAL(wp) ::  eta_s       !<
    7140     REAL(wp) ::  n           !< 3rd flattening
    7141     REAL(wp) ::  n2          !< n^2
    7142     REAL(wp) ::  n3          !< n^3
    7143     REAL(wp) ::  n4          !< n^4
    7144     REAL(wp) ::  n5          !< n^5
    7145     REAL(wp) ::  n6          !< n^6
    7146     REAL(wp) ::  nu          !<
    7147     REAL(wp) ::  nu_s        !<
    7148     REAL(wp) ::  sin_eta_s   !< sin(eta_s)
    7149     REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
    7150     REAL(wp) ::  tau_i       !<
    7151     REAL(wp) ::  tau_i_s     !<
    7152     REAL(wp) ::  tau_s       !<
    7153     REAL(wp) ::  x           !< adjusted easting
    7154     REAL(wp) ::  y           !< adjusted northing
    7155 
    7156     REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
    7157 
    7158     REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
    7159                                                !< (/semi_major_axis,
    7160                                                !<   inverse_flattening,
    7161                                                !<   longitude_of_prime_meridian,
    7162                                                !<   longitude_of_central_meridian,
    7163                                                !<   scale_factor_at_central_meridian,
    7164                                                !<   latitude_of_projection_origin,
    7165                                                !<   false_easting,
    7166                                                !<   false_northing /)
    7167 
    7168     x = eutm - crs(7)  ! remove false easting
    7169     y = nutm - crs(8)  ! remove false northing
    7170 
    7171 !-- from Karney 2011 Eq 15-22, 36:
    7172     e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
    7173     n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
    7174     n2 = n * n
    7175     n3 = n * n2
    7176     n4 = n * n3
    7177     n5 = n * n4
    7178     n6 = n * n5
    7179 
    7180     a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
    7181                                            + 0.015625_wp * n4   &
    7182                                            + 3.90625E-3_wp * n6 )
    7183 
    7184     nu  = x / ( crs(5) * a )
    7185     eta = y / ( crs(5) * a )
    7186 
    7187 !-- According to KrÃŒger (1912), eq. 26*
    7188     beta(1) =        0.5_wp                  * n  &
    7189             -        2.0_wp /         3.0_wp * n2 &
    7190             +       37.0_wp /        96.0_wp * n3 &
    7191             -        1.0_wp /       360.0_wp * n4 &
    7192             -       81.0_wp /       512.0_wp * n5 &
    7193             +    96199.0_wp /    604800.0_wp * n6
    7194 
    7195     beta(2) =        1.0_wp /        48.0_wp * n2 &
    7196             +        1.0_wp /        15.0_wp * n3 &
    7197             -      437.0_wp /      1440.0_wp * n4 &
    7198             +       46.0_wp /       105.0_wp * n5 &
    7199             -  1118711.0_wp /   3870720.0_wp * n6
    7200 
    7201     beta(3) =       17.0_wp /       480.0_wp * n3 &
    7202             -       37.0_wp /       840.0_wp * n4 &
    7203             -      209.0_wp /      4480.0_wp * n5 &
    7204             +     5569.0_wp /     90720.0_wp * n6
    7205 
    7206     beta(4) =     4397.0_wp /    161280.0_wp * n4 &
    7207             -       11.0_wp /       504.0_wp * n5 &
    7208             -   830251.0_wp /   7257600.0_wp * n6
    7209 
    7210     beta(5) =     4583.0_wp /    161280.0_wp * n5 &
    7211             -   108847.0_wp /   3991680.0_wp * n6
    7212 
    7213     beta(6) = 20648693.0_wp / 638668800.0_wp * n6
    7214 
    7215     eta_s = eta
    7216     nu_s  = nu
    7217     DO  j = 1, 6
    7218       eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
    7219       nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
    7220     ENDDO
    7221 
    7222     sinh_nu_s = SINH( nu_s )
    7223     sin_eta_s = SIN( eta_s )
    7224     cos_eta_s = COS( eta_s )
    7225 
    7226     tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
    7227 
    7228     tau_i = tau_s
    7229     delta_tau_i = 1.0_wp
    7230 
    7231     DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
    7232 
    7233       delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
    7234 
    7235       tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
    7236                - delta_i * SQRT( 1.0_wp + tau_i**2 )
    7237 
    7238       delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
    7239                    * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
    7240                    / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
    7241 
    7242       tau_i = tau_i + delta_tau_i
    7243 
    7244     ENDDO
    7245 
    7246     lat = ATAN( tau_i ) / pi * 180.0_wp
    7247     lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
    7248 
    7249  END SUBROUTINE convert_utm_to_geographic
    7250 
    72517106 END MODULE netcdf_interface
Note: See TracChangeset for help on using the changeset viewer.