| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121 |
- module m_confmap
- implicit none
- logical :: confmap_initialised = .false.
- real :: pi_1
- real :: pi_2
- real :: deg
- real :: rad
- real :: theta_a
- real :: phi_a
- real :: theta_b
- real :: phi_b
- real :: di
- real :: dj
- complex :: imagone
- complex :: ac
- complex :: bc
- complex :: cmna
- complex :: cmnb
- real :: mu_s
- real :: psi_s
- real :: epsil
- logical :: mercator
- real :: lat_a, lon_a
- real :: lat_b, lon_b
- real :: wlim, elim
- real :: slim, nlim
- real :: mercfac
- integer :: ires, jres
- contains
- ! This routine initializes constants used in the conformal mapping
- ! and must be called before the routines 'oldtonew' and 'newtoold'
- ! are called. The arguments of this routine are the locations of
- ! the two poles in the old coordiante system.
- !
- subroutine confmap_init(nx, ny)
- integer, intent(in) :: nx, ny
- real :: cx, cy, cz, theta_c, phi_c
- complex :: c, w
- logical :: ass, lold
-
- ! Read info file
- open(unit = 10, file = 'grid.info', form = 'formatted')
- read(10, *) lat_a, lon_a
- read(10, *) lat_b,lon_b
- read(10, *) wlim, elim, ires
- read(10, *) slim, nlim, jres
- read(10, *) ass
- read(10, *) ass
- read(10, *) ass
- read(10, *) mercator
- read(10, *) mercfac, lold
- close(10)
- if (ires /= nx .and. jres /= ny) then
- print *, 'initconfmap: WARNING -- the dimensions in grid.info are not'
- print *, 'initconfmap: WARNING -- consistent with nx and ny'
- print *, 'initconfmap: WARNING -- IGNORE IF RUNNING CURVIINT'
- stop '(initconfmap)'
- endif
- ! some constants
- !
- pi_1 = 3.14159265358979323846
- pi_2 = 0.5 * pi_1
- deg = 180.0 / pi_1
- rad = 1.0 / deg
- epsil = 1.0d-9
- di = (elim - wlim) / real(ires - 1) ! delta lon'
- dj = (nlim - slim) / real(jres - 1) ! delta lat' for spherical grid
- if (mercator) then
- dj = di
- if (lold) then
- print *, 'initconfmap: lold'
- slim = -mercfac * jres * dj
- else
- print *, 'initconfmap: not lold'
- slim = mercfac
- endif
- endif
- ! transform to spherical coordinates
- !
- theta_a = lon_a * rad
- phi_a = pi_2 - lat_a * rad
- theta_b = lon_b * rad
- phi_b = pi_2 - lat_b * rad
- ! find the angles of a vector pointing at a point located exactly
- ! between the poles
- !
- cx = cos(theta_a) * sin(phi_a) + cos(theta_b) * sin(phi_b)
- cy = sin(theta_a) * sin(phi_a) + sin(theta_b) * sin(phi_b)
- cz = cos(phi_a) + cos(phi_b)
- theta_c = atan2(cy, cx)
- phi_c = pi_2 - atan2(cz, sqrt(cx * cx + cy * cy))
- ! initialize constants used in the conformal mapping
- !
- imagone = (0.0, 1.0)
- ac = tan(0.5 * phi_a) * exp(imagone * theta_a)
- bc = tan(0.5 * phi_b) * exp(imagone * theta_b)
- c = tan(0.5 * phi_c) * exp(imagone * theta_c)
- cmna = c - ac
- cmnb = c - bc
- w = cmnb / cmna
- mu_s = atan2(aimag(w), real(w))
- psi_s = 2.0 * atan(abs(w))
- confmap_initialised = .true.
- end subroutine confmap_init
- end module m_confmap
|
