Changeset 4400 for palm/trunk/SOURCE/netcdf_interface_mod.f90
- Timestamp:
- Feb 10, 2020 8:32:41 PM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/netcdf_interface_mod.f90
r4360 r4400 25 25 ! ----------------- 26 26 ! $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 27 31 ! Adjusted output of multi-agent system for biometeorology 28 32 ! … … 135 139 136 140 USE netcdf_data_input_mod, & 137 ONLY: coord_ref_sys, init_model 141 ONLY: coord_ref_sys, & 142 crs_list, & 143 init_model 138 144 139 145 PRIVATE … … 417 423 418 424 USE basic_constants_and_equations_mod, & 419 ONLY: pi 425 ONLY: convert_utm_to_geographic, & 426 pi 420 427 421 428 USE control_parameters, & … … 564 571 565 572 REAL(wp), DIMENSION(1) :: last_time_coordinate !< last time value in file 566 REAL(wp), DIMENSION(8) :: crs_list !< list of coord_ref_sys values567 573 568 574 REAL(wp), DIMENSION(:), ALLOCATABLE :: netcdf_data !< … … 662 668 663 669 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 /)674 670 675 671 ! … … 7108 7104 7109 7105 7110 !------------------------------------------------------------------------------!7111 ! Description:7112 ! ------------7113 !> Convert UTM coordinates into geographic latitude and longitude. Conversion7114 !> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger287115 !> and Karney (2013) DOI: 10.1007/s00190-012-0578-z7116 !> Based on a JavaScript of the geodesy function library written by chrisveness7117 !> https://github.com/chrisveness/geodesy7118 !------------------------------------------------------------------------------!7119 SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )7120 7121 USE basic_constants_and_equations_mod, &7122 ONLY: pi7123 7124 IMPLICIT NONE7125 7126 INTEGER(iwp) :: j !< loop index7127 7128 REAL(wp), INTENT(in) :: eutm !< easting (UTM)7129 REAL(wp), INTENT(out) :: lat !< geographic latitude in degree7130 REAL(wp), INTENT(out) :: lon !< geographic longitude in degree7131 REAL(wp), INTENT(in) :: nutm !< northing (UTM)7132 7133 REAL(wp) :: a !< 2*pi*a is the circumference of a meridian7134 REAL(wp) :: cos_eta_s !< cos(eta_s)7135 REAL(wp) :: delta_i !<7136 REAL(wp) :: delta_tau_i !<7137 REAL(wp) :: e !< eccentricity7138 REAL(wp) :: eta !<7139 REAL(wp) :: eta_s !<7140 REAL(wp) :: n !< 3rd flattening7141 REAL(wp) :: n2 !< n^27142 REAL(wp) :: n3 !< n^37143 REAL(wp) :: n4 !< n^47144 REAL(wp) :: n5 !< n^57145 REAL(wp) :: n6 !< n^67146 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 easting7154 REAL(wp) :: y !< adjusted northing7155 7156 REAL(wp), DIMENSION(6) :: beta !< 6th order KrÃŒger expressions7157 7158 REAL(wp), DIMENSION(8), INTENT(in) :: crs !< coordinate reference system, consists of7159 !< (/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 easting7169 y = nutm - crs(8) ! remove false northing7170 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 * n7175 n3 = n * n27176 n4 = n * n37177 n5 = n * n47178 n6 = n * n57179 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 * n67194 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 * n67200 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 * n67205 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 * n67209 7210 beta(5) = 4583.0_wp / 161280.0_wp * n5 &7211 - 108847.0_wp / 3991680.0_wp * n67212 7213 beta(6) = 20648693.0_wp / 638668800.0_wp * n67214 7215 eta_s = eta7216 nu_s = nu7217 DO j = 1, 67218 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 ENDDO7221 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_s7229 delta_tau_i = 1.0_wp7230 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_i7243 7244 ENDDO7245 7246 lat = ATAN( tau_i ) / pi * 180.0_wp7247 lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)7248 7249 END SUBROUTINE convert_utm_to_geographic7250 7251 7106 END MODULE netcdf_interface
Note: See TracChangeset
for help on using the changeset viewer.