!***************************************************************************************** !> ! IGRF model ! !### History ! * SHELLIG.FOR, Version 2.0, January 1992 ! * 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2 ! * 1/27/92-DKB- Adopted to IGRF-91 coefficients model ! * 2/05/92-DKB- Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE ! * 8/08/95-DKB- Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S ! * 5/31/00-DKB- Updated to IGRF-45-00; new coeff.: IGRF00, IGRF00s ! * 3/24/05-DKB- Updated to IGRF-45-10; new coeff.: IGRF05, IGRF05s module shellig_module use radbelt_kinds_module implicit none private integer, parameter :: filename_len = 14 !! length of the model data file names ! parameters formerly in `gener` common block real(wp), parameter :: Era = 6371.2_wp !! earth radius for normalization of cartesian coordinates (6371.2 km) real(wp), parameter :: erequ = 6378.16_wp real(wp), parameter :: erpol = 6356.775_wp real(wp), parameter :: Aquad = erequ * erequ !! square of major half axis for !! earth ellipsoid as recommended by international !! astronomical union real(wp), parameter :: Bquad = erpol * erpol !! square of minor half axis for !! earth ellipsoid as recommended by international !! astronomical union real(wp), parameter :: Umr = atan(1.0_wp) * 4.0_wp / 180.0_wp !! atan(1.0)*4./180. <degree>*umr=<radiant> real(wp), dimension(3, 3), parameter :: u = reshape([+0.3511737_wp, -0.9148385_wp, -0.1993679_wp, & +0.9335804_wp, +0.3583680_wp, +0.0000000_wp, & +0.0714471_wp, -0.1861260_wp, +0.9799247_wp], [3, 3]) integer, parameter :: max_loop_index = 3333 !! used in [[shellg]] for the field line integration loop type, public :: shellig_type private character(len=:), allocatable :: igrf_dir !! directory containing the data files ! formerly in the `fidb0` common block real(wp), dimension(3) :: sp = 0.0_wp ! formerly in blank common real(wp), dimension(3) :: xi = 0.0_wp real(wp), dimension(144) :: h = 0.0_wp !! Field model coefficients adjusted for [[shellg]] ! formerly in `model` common block integer :: iyea = 0 !! the int year corresponding to the file `name` that has been read character(len=:), allocatable :: name !! file name integer :: nmax = 0 !! maximum order of spherical harmonics real(wp) :: Time = 0.0_wp !! year (decimal: 1973.5) for which magnetic field is to be calculated real(wp), dimension(144) :: g = 0.0_wp !! `g(m)` -- normalized field coefficients (see [[feldcof]]) m=nmax*(nmax+2) integer :: nmax1 = 0 !! saved variables from the file integer :: nmax2 = 0 !! saved variables from the file real(wp), dimension(144) :: g_cache = 0.0_wp !! saved `g` from the file ! formerly saved vars in shellg: real(wp) :: step = 0.20_wp !! step size for field line tracing real(wp) :: steq = 0.03_wp !! step size for integration ! from feldcof, so we can cache the coefficients real(wp), dimension(120) :: gh2 = 0.0_wp ! JW : why is this 120 and g is 144 ??? real(wp), dimension(:, :), allocatable :: p !! this was `p(8,100)` in the original code. !! used for the field line integration loop. !! changed it to be allocatable since it was !! changed to be p(8,3334). contains private procedure, public :: igrf, igrfc procedure, public :: feldcof procedure, public :: feldg, feldc procedure, public :: shellg, shellc procedure, public :: findb0 procedure :: stoer, feldi procedure, public :: set_data_file_dir, get_data_file_dir procedure, public :: destroy => destroy_shellig_type end type shellig_type contains !***************************************************************************************** !***************************************************************************************** !> ! Destroy a [[shellig_type]]. subroutine destroy_shellig_type(me) class(shellig_type), intent(out) :: me end subroutine destroy_shellig_type !***************************************************************************************** !> ! Set the directory containing the data files. subroutine set_data_file_dir(me, dir) class(shellig_type), intent(inout) :: me character(len=*), intent(in) :: dir me%igrf_dir = trim(dir) end subroutine set_data_file_dir !***************************************************************************************** !> ! Get the directory containing the data files. function get_data_file_dir(me) result(dir) class(shellig_type), intent(in) :: me character(len=:), allocatable :: dir if (allocated(me%igrf_dir)) then dir = trim(me%igrf_dir)//'/' else dir = 'data/igrf/' ! default end if end function get_data_file_dir !***************************************************************************************** !> ! Wrapper for IGRF functions. subroutine igrf(me, lon, lat, height, year, xl, bbx) class(shellig_type), intent(inout) :: me real(wp), intent(in) :: lon !! geodetic longitude in degrees (east) real(wp), intent(in) :: lat !! geodetic latitude in degrees (north) real(wp), intent(in) :: height !! altitude in km above sea level real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) real(wp), intent(out) :: xl !! l-value real(wp), intent(out) :: bbx !! b_total / b_equatorial ratio real(wp) :: bab1, babs, bdel, bdown, beast, & beq, bequ, bnorth, dimo, rr0 integer :: icode logical :: val real(wp), parameter :: stps = 0.05_wp ! JW : do we need to reset some or all of these ? me%sp = 0.0_wp me%xi = 0.0_wp me%h = 0.0_wp me%step = 0.20_wp me%steq = 0.03_wp call me%feldcof(year, dimo) call me%feldg(lat, lon, height, bnorth, beast, bdown, babs) call me%shellg(lat, lon, height, dimo, xl, icode, bab1) bequ = dimo / (xl * xl * xl) if (icode == 1) then bdel = 1.0e-3_wp call me%findb0(stps, bdel, val, beq, rr0) if (val) bequ = beq end if bbx = babs / bequ end subroutine igrf !***************************************************************************************** !***************************************************************************************** !> ! Alternate version of [[igrf]] for cartesian coordinates. subroutine igrfc(me, v, year, xl, bbx) class(shellig_type), intent(inout) :: me real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! x-axis pointing to equator at 0 longitude !! y-axis pointing to equator at 90 long. !! z-axis pointing to north pole real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) real(wp), intent(out) :: xl !! l-value real(wp), intent(out) :: bbx !! b_total / b_equatorial ratio real(wp) :: bab1, bdel, beq, bequ, dimo, rr0 integer :: icode logical :: val real(wp), dimension(3) :: b real(wp), parameter :: stps = 0.05_wp ! JW : do we need to reset some or all of these ? me%sp = 0.0_wp me%xi = 0.0_wp me%h = 0.0_wp me%step = 0.20_wp me%steq = 0.03_wp call me%feldcof(year, dimo) call me%feldc(v, b) call me%shellc(v, dimo, xl, icode, bab1) bequ = dimo / (xl * xl * xl) if (icode == 1) then bdel = 1.0e-3_wp call me%findb0(stps, bdel, val, beq, rr0) if (val) bequ = beq end if bbx = norm2(b) / bequ end subroutine igrfc !***************************************************************************************** !***************************************************************************************** !> subroutine findb0(me, stps, bdel, value, bequ, rr0) class(shellig_type), intent(inout) :: me real(wp), intent(in) :: stps real(wp), intent(inout) :: bdel real(wp), intent(out) :: bequ logical, intent(out) :: value real(wp), intent(out) :: rr0 real(wp) :: b, bdelta, bmin, bold, bq1, & bq2, bq3, p(8, 4), r1, r2, r3, & rold, step, step12, zz integer :: i, irun, j, n step = stps irun = 0 rold = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings main: do irun = irun + 1 if (irun > 5) then value = .false. exit main end if ! first three points p(1, 2) = me%sp(1) p(2, 2) = me%sp(2) p(3, 2) = me%sp(3) step = -sign(step, p(3, 2)) call me%stoer(p(1, 2), bq2, r2) p(1, 3) = p(1, 2) + 0.5_wp * step * p(4, 2) p(2, 3) = p(2, 2) + 0.5_wp * step * p(5, 2) p(3, 3) = p(3, 2) + 0.5_wp * step call me%stoer(p(1, 3), bq3, r3) p(1, 1) = p(1, 2) - step * (2.0_wp * p(4, 2) - p(4, 3)) p(2, 1) = p(2, 2) - step * (2.0_wp * p(5, 2) - p(5, 3)) p(3, 1) = p(3, 2) - step call me%stoer(p(1, 1), bq1, r1) p(1, 3) = p(1, 2) + step * (20.0_wp * p(4, 3) - 3.*p(4, 2) + p(4, 1)) / 18.0_wp p(2, 3) = p(2, 2) + step * (20.0_wp * p(5, 3) - 3.*p(5, 2) + p(5, 1)) / 18.0_wp p(3, 3) = p(3, 2) + step call me%stoer(p(1, 3), bq3, r3) ! invert sense if required if (bq3 > bq1) then step = -step r3 = r1 bq3 = bq1 do i = 1, 5 zz = p(i, 1) p(i, 1) = p(i, 3) p(i, 3) = zz end do end if ! initialization step12 = step / 12.0_wp value = .true. bmin = 1.0e4_wp bold = 1.0e4_wp ! corrector (field line tracing) n = 0 corrector: do p(1, 3) = p(1, 2) + step12 * (5.0_wp * p(4, 3) + 8.0_wp * p(4, 2) - p(4, 1)) n = n + 1 p(2, 3) = p(2, 2) + step12 * (5.0_wp * p(5, 3) + 8.0_wp * p(5, 2) - p(5, 1)) ! predictor (field line tracing) p(1, 4) = p(1, 3) + step12 * (23.0_wp * p(4, 3) - 16.0_wp * p(4, 2) + 5.0_wp * p(4, 1)) p(2, 4) = p(2, 3) + step12 * (23.0_wp * p(5, 3) - 16.0_wp * p(5, 2) + 5.0_wp * p(5, 1)) p(3, 4) = p(3, 3) + step call me%stoer(p(1, 4), bq3, r3) do j = 1, 3 do i = 1, 8 p(i, j) = p(i, j + 1) end do end do b = sqrt(bq3) if (b < bmin) bmin = b if (b > bold) exit corrector bold = b rold = 1.0_wp / r3 me%sp(1) = p(1, 4) me%sp(2) = p(2, 4) me%sp(3) = p(3, 4) end do corrector if (bold /= bmin) value = .false. bdelta = (b - bold) / bold if (bdelta <= bdel) exit main step = step / 10.0_wp end do main rr0 = rold bequ = bold bdel = bdelta end subroutine findb0 !***************************************************************************************** !> ! Wrapper to [[shellg]] to be used with cartesian coordinates. ! !@note In the original code, this was an ENTRY point in [[shellg]] and didn't ! include all the outputs. subroutine shellc(me, v, dimo, fl, icode, b0) class(shellig_type), intent(inout) :: me real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. !! * z-axis pointing to north pole real(wp), intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) real(wp), intent(out) :: fl !! l-value integer, intent(out) :: icode !! * =1 normal completion !! * =2 unphysical conjugate point (fl meaningless) !! * =3 shell parameter greater than limit up to !! which accurate calculation is required; !! approximation is used. real(wp), intent(out) :: b0 !! magnetic field strength in gauss real(wp) :: glat, glon, alt !! not used call me%shellg(glat, glon, alt, dimo, fl, icode, b0, v) end subroutine shellc !***************************************************************************************** !> ! calculates l-value for specified geodaetic coordinates, altitude ! and gemagnetic field model. ! !### Reference ! * G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE ! NO. 67, 1970. ! * G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972 ! !### History ! * CHANGES (D. BILITZA, NOV 87): ! - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/ ! - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990 subroutine shellg(me, glat, glon, alt, dimo, fl, icode, b0, v) class(shellig_type), intent(inout) :: me real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) real(wp), intent(in) :: alt !! altitude in km above sea level real(wp), intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) real(wp), intent(out) :: fl !! l-value integer, intent(out) :: icode !! * =1 normal completion !! * =2 unphysical conjugate point (fl meaningless) !! * =3 shell parameter greater than limit up to !! which accurate calculation is required; !! approximation is used. real(wp), intent(out) :: b0 !! magnetic field strength in gauss real(wp), dimension(3), intent(in), optional :: v !! cartesian coordinates in earth radii (6371.2 km) !! !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. !! * z-axis pointing to north pole !! !! If this argument is present, it is used !! instead of glat,glon,alt. See [[shellc]]. real(wp) :: arg1, arg2, bequ, bq1, bq2, bq3, c0, c1, c2, c3, & d0, d1, d2, dimob0, e0, e1, e2, ff, fi, gg, & hli, oradik, oterm, r, r1, r2, r3, r3h, radik, & rq, step12, step2, stp, t, term, xx, z, zq, zz integer :: i, iequ, n real(wp), parameter :: rmin = 0.05_wp !! boundaries for identification of `icode=2 and 3` real(wp), parameter :: rmax = 1.01_wp !! boundaries for identification of `icode=2 and 3` if (.not. allocated(me%p)) allocate (me%p(8, max_loop_index + 1)) ! because `p(:,n+1)` in the loop bequ = 1.0e10_wp if (present(v)) then me%xi(1) = v(1) me%xi(2) = v(2) me%xi(3) = v(3) else me%xi = geo_to_cart(glat, glon, alt) end if associate (p => me%p) radik = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings ! convert to dipol-oriented co-ordinates rq = 1.0_wp / (me%xi(1) * me%xi(1) + me%xi(2) * me%xi(2) + me%xi(3) * me%xi(3)) r3h = sqrt(rq * sqrt(rq)) p(1, 2) = (me%xi(1) * u(1, 1) + me%xi(2) * u(2, 1) + me%xi(3) * u(3, 1)) * r3h p(2, 2) = (me%xi(1) * u(1, 2) + me%xi(2) * u(2, 2)) * r3h p(3, 2) = (me%xi(1) * u(1, 3) + me%xi(2) * u(2, 3) + me%xi(3) * u(3, 3)) * rq ! first three points of field line me%step = -sign(me%step, p(3, 2)) call me%stoer(p(1, 2), bq2, r2) b0 = sqrt(bq2) p(1, 3) = p(1, 2) + 0.5_wp * me%step * p(4, 2) p(2, 3) = p(2, 2) + 0.5_wp * me%step * p(5, 2) p(3, 3) = p(3, 2) + 0.5_wp * me%step call me%stoer(p(1, 3), bq3, r3) p(1, 1) = p(1, 2) - me%step * (2.0_wp * p(4, 2) - p(4, 3)) p(2, 1) = p(2, 2) - me%step * (2.0_wp * p(5, 2) - p(5, 3)) p(3, 1) = p(3, 2) - me%step call me%stoer(p(1, 1), bq1, r1) p(1, 3) = p(1, 2) + me%step * (20.0_wp * p(4, 3) - 3.*p(4, 2) + p(4, 1)) / 18.0_wp p(2, 3) = p(2, 2) + me%step * (20.0_wp * p(5, 3) - 3.*p(5, 2) + p(5, 1)) / 18.0_wp p(3, 3) = p(3, 2) + me%step call me%stoer(p(1, 3), bq3, r3) ! invert sense if required if (bq3 > bq1) then me%step = -me%step r3 = r1 bq3 = bq1 do i = 1, 7 zz = p(i, 1) p(i, 1) = p(i, 3) p(i, 3) = zz end do end if ! search for lowest magnetic field strength if (bq1 < bequ) then bequ = bq1 iequ = 1 end if if (bq2 < bequ) then bequ = bq2 iequ = 2 end if if (bq3 < bequ) then bequ = bq3 iequ = 3 end if ! initialization of integration loops step12 = me%step / 12.0_wp step2 = me%step + me%step me%steq = sign(me%steq, me%step) fi = 0.0_wp icode = 1 oradik = 0.0_wp oterm = 0.0_wp stp = r2 * me%steq z = p(3, 2) + stp stp = stp / 0.75_wp p(8, 1) = step2 * (p(1, 1) * p(4, 1) + p(2, 1) * p(5, 1)) p(8, 2) = step2 * (p(1, 2) * p(4, 2) + p(2, 2) * p(5, 2)) ! main loop (field line tracing) main: do n = 3, max_loop_index ! corrector (field line tracing) p(1, n) = p(1, n - 1) + step12 * (5.0_wp * p(4, n) + 8.0_wp * p(4, n - 1) - p(4, n - 2)) p(2, n) = p(2, n - 1) + step12 * (5.0_wp * p(5, n) + 8.0_wp * p(5, n - 1) - p(5, n - 2)) ! prepare expansion coefficients for interpolation ! of slowly varying quantities p(8, n) = step2 * (p(1, n) * p(4, n) + p(2, n) * p(5, n)) c0 = p(1, n - 1)**2 + p(2, n - 1)**2 c1 = p(8, n - 1) c2 = (p(8, n) - p(8, n - 2)) * 0.25_wp c3 = (p(8, n) + p(8, n - 2) - c1 - c1) / 6.0_wp d0 = p(6, n - 1) d1 = (p(6, n) - p(6, n - 2)) * 0.5_wp d2 = (p(6, n) + p(6, n - 2) - d0 - d0) * 0.5_wp e0 = p(7, n - 1) e1 = (p(7, n) - p(7, n - 2)) * 0.5_wp e2 = (p(7, n) + p(7, n - 2) - e0 - e0) * 0.5_wp inner: do ! inner loop (for quadrature) t = (z - p(3, n - 1)) / me%step if (t > 1.0_wp) then ! predictor (field line tracing) p(1, n + 1) = p(1, n) + step12 * (23.0_wp * p(4, n) - 16.0_wp * p(4, n - 1) + 5.0_wp * p(4, n - 2)) p(2, n + 1) = p(2, n) + step12 * (23.0_wp * p(5, n) - 16.0_wp * p(5, n - 1) + 5.0_wp * p(5, n - 2)) p(3, n + 1) = p(3, n) + me%step call me%stoer(p(1, n + 1), bq3, r3) ! search for lowest magnetic field strength if (bq3 < bequ) then iequ = n + 1 bequ = bq3 end if exit inner else hli = 0.5_wp * (((c3 * t + c2) * t + c1) * t + c0) zq = z * z r = hli + sqrt(hli * hli + zq) if (r <= rmin) then ! approximation for high values of l. icode = 3 t = -p(3, n - 1) / me%step fl = 1.0_wp / (abs(((c3 * t + c2) * t + c1) * t + c0) + 1.0e-15_wp) return end if rq = r * r ff = sqrt(1.0_wp + 3.0_wp * zq / rq) radik = b0 - ((d2 * t + d1) * t + d0) * r * rq * ff if (r > rmax) then icode = 2 radik = radik - 12.0_wp * (r - rmax)**2 end if if (radik + radik <= oradik) exit main term = sqrt(radik) * ff * ((e2 * t + e1) * t + e0) / (rq + zq) fi = fi + stp * (oterm + term) oradik = radik oterm = term stp = r * me%steq z = z + stp end if end do inner end do main if (iequ < 2) iequ = 2 me%sp(1) = p(1, iequ - 1) me%sp(2) = p(2, iequ - 1) me%sp(3) = p(3, iequ - 1) if (oradik >= 1.0e-15_wp) fi = fi + stp / & 0.75_wp * oterm * oradik / & (oradik - radik) ! the minimal allowable value of fi was changed from 1e-15 to 1e-12, ! because 1e-38 is the minimal allowable arg. for alog in our envir. ! d. bilitza, nov 87. fi = 0.5_wp * abs(fi) / sqrt(b0) + 1.0e-12_wp ! compute l from b and i. same as carmel in invar. ! correct dipole moment is used here. d. bilitza, nov 87. dimob0 = dimo / b0 arg1 = log(fi) arg2 = log(dimob0) ! arg = fi*fi*fi/dimob0 ! if(abs(arg)>88.0_wp) arg=88.0_wp xx = 3 * arg1 - arg2 if (xx > 23.0_wp) then gg = xx - 3.0460681_wp elseif (xx > 11.7_wp) then gg = (((((2.8212095e-8_wp * xx - 3.8049276e-6_wp) * xx + & 2.170224e-4_wp) * xx - 6.7310339e-3_wp) * xx + & 1.2038224e-1_wp) * xx - 1.8461796e-1_wp) * xx + 2.0007187_wp elseif (xx > +3.0_wp) then gg = ((((((((6.3271665e-10_wp * xx - 3.958306e-8_wp) * xx + & 9.9766148e-07_wp) * xx - 1.2531932e-5_wp) * xx + & 7.9451313e-5_wp) * xx - 3.2077032e-4_wp) * xx + & 2.1680398e-3_wp) * xx + 1.2817956e-2_wp) * xx + & 4.3510529e-1_wp) * xx + 6.222355e-1_wp elseif (xx > -3.0_wp) then gg = ((((((((2.6047023e-10_wp * xx + 2.3028767e-9_wp) * xx - & 2.1997983e-8_wp) * xx - 5.3977642e-7_wp) * xx - & 3.3408822e-6_wp) * xx + 3.8379917e-5_wp) * xx + & 1.1784234e-3_wp) * xx + 1.4492441e-2_wp) * xx + & 4.3352788e-1_wp) * xx + 6.228644e-1_wp elseif (xx > -22.0_wp) then gg = ((((((((-8.1537735e-14_wp * xx + 8.3232531e-13_wp) * xx + & 1.0066362e-9_wp) * xx + 8.1048663e-8_wp) * xx + & 3.2916354e-6_wp) * xx + 8.2711096e-5_wp) * xx + & 1.3714667e-3_wp) * xx + 1.5017245e-2_wp) * xx + & 4.3432642e-1_wp) * xx + 6.2337691e-1_wp else gg = 3.33338e-1_wp * xx + 3.0062102e-1_wp end if fl = exp(log((1.0_wp + exp(gg)) * dimob0) / 3.0_wp) end associate end subroutine shellg !***************************************************************************************** !> ! subroutine used for field line tracing in [[shellg]]. ! calls entry point [[feldi]] in geomagnetic field subroutine [[feldg]] subroutine stoer(me, p, bq, r) class(shellig_type), intent(inout) :: me real(wp), dimension(7), intent(inout) :: p real(wp), intent(out) :: bq real(wp), intent(out) :: r real(wp) :: dr, dsq, dx, dxm, dy, dym, dz, & dzm, fli, q, rq, wr, xm, ym, zm ! xm,ym,zm are geomagnetic cartesian inverse co-ordinates zm = P(3) fli = P(1) * P(1) + P(2) * P(2) + 1.0e-15_wp R = 0.5_wp * (fli + sqrt(fli * fli + (zm + zm)**2)) rq = R * R wr = sqrt(R) xm = P(1) * wr ym = P(2) * wr ! transform to geographic co-ordinate system me%Xi(1) = xm * u(1, 1) + ym * u(1, 2) + zm * u(1, 3) me%Xi(2) = xm * u(2, 1) + ym * u(2, 2) + zm * u(2, 3) me%Xi(3) = xm * u(3, 1) + zm * u(3, 3) ! compute derivatives ! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results ! are the same; dkb Feb 1998. ! JW : feb 2024 : xi, h now class variables. call me%feldi() q = me%H(1) / rq dx = me%H(3) + me%H(3) + q * me%Xi(1) dy = me%H(4) + me%H(4) + q * me%Xi(2) dz = me%H(2) + me%H(2) + q * me%Xi(3) ! transform back to geomagnetic co-ordinate system dxm = u(1, 1) * dx + u(2, 1) * dy + u(3, 1) * dz dym = u(1, 2) * dx + u(2, 2) * dy dzm = u(1, 3) * dx + u(2, 3) * dy + u(3, 3) * dz dr = (xm * dxm + ym * dym + zm * dzm) / R ! form slowly varying expressions P(4) = (wr * dxm - 0.5_wp * P(1) * dr) / (R * dzm) P(5) = (wr * dym - 0.5_wp * P(2) * dr) / (R * dzm) dsq = rq * (dxm * dxm + dym * dym + dzm * dzm) Bq = dsq * rq * rq P(6) = sqrt(dsq / (rq + 3.0_wp * zm * zm)) P(7) = P(6) * (rq + zm * zm) / (rq * dzm) end subroutine stoer !***************************************************************************************** !> ! Calculates earth magnetic field from spherical harmonics model ! !### Reference ! ref: g. kluge, european space operations centre, internal note 61, ! 1970. ! !### History ! * changes (d. bilitza, nov 87): ! - field coefficients in binary data files instead of block data ! - calculates dipol moment ! !@note In the original code, [[feldc] and [[feldi]] were ! ENTRY points to this routine subroutine feldg(me, glat, glon, alt, bnorth, beast, bdown, babs) class(shellig_type), intent(inout) :: me real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) real(wp), intent(in) :: alt !! altitude in km above sea level real(wp), intent(out) :: bnorth, beast, bdown !! components of the field with respect !! to the local geodetic coordinate system, with axis !! pointing in the tangential plane to the north, east !! and downward. real(wp), intent(out) :: Babs !! magnetic field strength in gauss real(wp) :: brho, bxxx, byyy, bzzz, cp, ct, d, f, rho, & rlat, rlon, rq, s, sp, st, t, & x, xxx, y, yyy, z, zzz integer :: i, ih, ihmax, il, imax, k, last, m ! same calculation as geo_to_cart, but not used here ! because the intermediate variables are also used below. rlat = glat * umr ct = sin(rlat) st = cos(rlat) d = sqrt(aquad - (aquad - bquad) * ct * ct) rlon = glon * umr cp = cos(rlon) sp = sin(rlon) zzz = (alt + bquad / d) * ct / era rho = (alt + aquad / d) * st / era xxx = rho * cp yyy = rho * sp rq = 1.0_wp / (xxx * xxx + yyy * yyy + zzz * zzz) me%xi = [xxx, yyy, zzz] * rq ihmax = me%nmax * me%nmax + 1 last = ihmax + me%nmax + me%nmax imax = me%nmax + me%nmax - 1 do i = ihmax, last me%h(i) = me%g(i) end do do k = 1, 3, 2 i = imax ih = ihmax do il = ih - i f = 2.0_wp / real(i - k + 2, wp) x = me%xi(1) * f y = me%xi(2) * f z = me%xi(3) * (f + f) i = i - 2 if ((i - 1) >= 0) then if ((i - 1) > 0) then do m = 3, i, 2 me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) end do end if me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) end if me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) ih = il if (i < k) exit end do end do s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2)) t = (rq + rq) * sqrt(rq) bxxx = t * (me%h(3) - s * xxx) byyy = t * (me%h(4) - s * yyy) bzzz = t * (me%h(2) - s * zzz) babs = sqrt(bxxx * bxxx + byyy * byyy + bzzz * bzzz) beast = byyy * cp - bxxx * sp brho = byyy * sp + bxxx * cp bnorth = bzzz * st - brho * ct bdown = -bzzz * ct - brho * st end subroutine feldg !***************************************************************************************** !> ! Alternate version of [[feldg]] to be used with cartesian coordinates subroutine feldc(me, v, b) class(shellig_type), intent(inout) :: me real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) !! x-axis pointing to equator at 0 longitude !! y-axis pointing to equator at 90 long. !! z-axis pointing to north pole real(wp), intent(out) :: b(3) !! field components real(wp) :: f, rq, s, t, x, xxx, y, yyy, z, zzz integer :: i, ih, ihmax, il, imax, k, last, m xxx = v(1) yyy = v(2) zzz = v(3) rq = 1.0_wp / (xxx * xxx + yyy * yyy + zzz * zzz) me%xi = [xxx, yyy, zzz] * rq ihmax = me%nmax * me%nmax + 1 last = ihmax + me%nmax + me%nmax imax = me%nmax + me%nmax - 1 do i = ihmax, last me%h(i) = me%g(i) end do do k = 1, 3, 2 i = imax ih = ihmax do il = ih - i f = 2.0_wp / real(i - k + 2, wp) x = me%xi(1) * f y = me%xi(2) * f z = me%xi(3) * (f + f) i = i - 2 if ((i - 1) >= 0) then if ((i - 1) > 0) then do m = 3, i, 2 me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) end do end if me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) end if me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) ih = il if (i < k) exit end do end do s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2)) t = (rq + rq) * sqrt(rq) b(1) = t * (me%h(3) - s * xxx) b(2) = t * (me%h(4) - s * yyy) b(3) = t * (me%h(2) - s * zzz) end subroutine feldc !***************************************************************************************** !> ! Used for `l` computation. subroutine feldi(me) class(shellig_type), intent(inout) :: me real(wp) :: f, x, y, z integer :: i, ih, ihmax, il, imax, k, last, m ihmax = me%nmax * me%nmax + 1 last = ihmax + me%nmax + me%nmax imax = me%nmax + me%nmax - 1 do i = ihmax, last me%h(i) = me%g(i) end do do k = 1, 3, 2 i = imax ih = ihmax do il = ih - i f = 2.0_wp / real(i - k + 2, wp) x = me%xi(1) * f y = me%xi(2) * f z = me%xi(3) * (f + f) i = i - 2 if ((i - 1) >= 0) then if ((i - 1) > 0) then do m = 3, i, 2 me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - & me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2)) me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - & me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1)) end do end if me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih)) me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih)) end if me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2)) ih = il if (i < k) exit end do end do end subroutine feldi !***************************************************************************************** !> ! Determines coefficients and dipol moment from IGRF models ! !### Author ! * D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771, ! (301) 286-9536 NOV 1987. ! !### History ! * corrected for 2000 update - dkb- 5/31/2000 ! * updated to IGRF-2000 version -dkb- 5/31/2000 ! * updated to IGRF-2005 version -dkb- 3/24/2000 subroutine feldcof(me, year, dimo) class(shellig_type), intent(inout) :: me real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to !! be calculated (e.g.:1995.5 for day 185 of 1995) real(wp), intent(out) :: dimo !! geomagnetic dipol moment in gauss (normalized !! to earth's radius) at the time (year) real(wp) :: dte1, dte2, erad, gha(144), sqrt2 integer :: i, ier, j, l, m, n, iyea character(len=:), allocatable :: fil2 real(wp) :: x, f0, f !! these were double precision in original !! code while everything else was single precision ! changed to conform with IGRF 45-95, also FILMOD, DTEMOD arrays +1 character(len=filename_len), dimension(17), parameter :: filmod = [ & 'dgrf1945.dat ', 'dgrf1950.dat ', 'dgrf1955.dat ', 'dgrf1960.dat ', & 'dgrf1965.dat ', 'dgrf1970.dat ', 'dgrf1975.dat ', 'dgrf1980.dat ', & 'dgrf1985.dat ', 'dgrf1990.dat ', 'dgrf1995.dat ', 'dgrf2000.dat ', & 'dgrf2005.dat ', 'dgrf2010.dat ', 'dgrf2015.dat ', 'igrf2020.dat ', & 'igrf2020s.dat'] real(wp), dimension(17), parameter :: dtemod = [1945.0_wp, 1950.0_wp, 1955.0_wp, & 1960.0_wp, 1965.0_wp, 1970.0_wp, & 1975.0_wp, 1980.0_wp, 1985.0_wp, & 1990.0_wp, 1995.0_wp, 2000.0_wp, & 2005.0_wp, 2010.0_wp, 2015.0_wp, & 2020.0_wp, 2025.0_wp] integer, parameter :: numye = size(dtemod) - 1 ! number of 5-year priods represented by IGRF integer, parameter :: is = 0 !! * is=0 for schmidt normalization !! * is=1 gauss normalization logical :: read_file !-- determine igrf-years for input-year me%time = year iyea = int(year / 5.0_wp) * 5 read_file = iyea /= me%iyea ! if we have to read the file me%iyea = iyea l = (me%iyea - 1945) / 5 + 1 if (l < 1) l = 1 if (l > numye) l = numye dte1 = dtemod(l) me%name = me%get_data_file_dir()//trim(filmod(l)) dte2 = dtemod(l + 1) fil2 = me%get_data_file_dir()//trim(filmod(l + 1)) if (read_file) then ! get igrf coefficients for the boundary years ! [if they have not ready been loaded] call getshc(me%name, me%nmax1, erad, me%g, ier) if (ier /= 0) error stop 'error reading file: '//trim(me%name) me%g_cache = me%g ! because it is modified below, we have to cache the original values from the file call getshc(fil2, me%nmax2, erad, me%gh2, ier) if (ier /= 0) error stop 'error reading file: '//trim(fil2) else me%g = me%g_cache end if !-- determine igrf coefficients for year if (l <= numye - 1) then call intershc(year, dte1, me%nmax1, me%g, dte2, me%nmax2, me%gh2, me%nmax, gha) else call extrashc(year, dte1, me%nmax1, me%g, me%nmax2, me%gh2, me%nmax, gha) end if !-- determine magnetic dipol moment and coeffiecients g f0 = 0.0_wp do j = 1, 3 f = gha(j) * 1.0e-5_wp f0 = f0 + f * f end do dimo = sqrt(f0) me%g(1) = 0.0_wp i = 2 f0 = 1.0e-5_wp if (is == 0) f0 = -f0 sqrt2 = sqrt(2.0_wp) do n = 1, me%nmax x = n f0 = f0 * x * x / (4.0_wp * x - 2.0_wp) if (is == 0) f0 = f0 * (2.0_wp * x - 1.0_wp) / x f = f0 * 0.5_wp if (is == 0) f = f * sqrt2 me%g(i) = gha(i - 1) * f0 i = i + 1 do m = 1, n f = f * (x + m) / (x - m + 1.0_wp) if (is == 0) f = f * sqrt((x - m + 1.0_wp) / (x + m)) me%g(i) = gha(i - 1) * f me%g(i + 1) = gha(i) * f i = i + 2 end do end do end subroutine feldcof !***************************************************************************************** !> ! Reads spherical harmonic coefficients from the specified ! file into an array. ! !### Author ! * Version 1.01, A. Zunde, USGS, MS 964, ! Box 25046 Federal Center, Denver, CO 80225 subroutine getshc(Fspec, Nmax, Erad, Gh, Ier) character(len=*), intent(in) :: Fspec !! File specification integer, intent(out) :: Nmax !! Maximum degree and order of model real(wp), intent(out) :: Erad !! Earth's radius associated with the spherical !! harmonic coefficients, in the same units as !! elevation real(wp), dimension(*), intent(out) :: Gh !! Schmidt quasi-normal internal spherical !! harmonic coefficients integer, intent(out) :: Ier !! Error number: !! !! * 0, no error !! * -2, records out of order !! * FORTRAN run-time error number integer :: iu !! logical unit number real(wp) :: g, h integer :: i, m, mm, n, nn read_file: block ! --------------------------------------------------------------- ! Open coefficient file. Read past first header record. ! Read degree and order of model and Earth's radius. ! --------------------------------------------------------------- open (newunit=Iu, FILE=Fspec, STATUS='OLD', IOSTAT=Ier) if (Ier /= 0) then write (*, *) 'Error opening file: '//trim(fspec) exit read_file end if read (Iu, *, IOSTAT=Ier) if (Ier /= 0) exit read_file read (Iu, *, IOSTAT=Ier) Nmax, Erad if (Ier /= 0) exit read_file ! --------------------------------------------------------------- ! Read the coefficient file, arranged as follows: ! ! N M G H ! ---------------------- ! / 1 0 GH(1) - ! / 1 1 GH(2) GH(3) ! / 2 0 GH(4) - ! / 2 1 GH(5) GH(6) ! NMAX*(NMAX+3)/2 / 2 2 GH(7) GH(8) ! records \ 3 0 GH(9) - ! \ . . . . ! \ . . . . ! NMAX*(NMAX+2) \ . . . . ! elements in GH \ NMAX NMAX . . ! ! N and M are, respectively, the degree and order of the ! coefficient. ! --------------------------------------------------------------- i = 0 main: do nn = 1, Nmax do mm = 0, nn read (Iu, *, IOSTAT=Ier) n, m, g, h if (Ier /= 0) exit main if (nn /= n .or. mm /= m) then Ier = -2 exit main end if i = i + 1 Gh(i) = g if (m /= 0) then i = i + 1 Gh(i) = h end if end do end do main end block read_file close (Iu) end subroutine getshc !***************************************************************************************** !> ! Interpolates linearly, in time, between two spherical ! harmonic models. ! ! The coefficients (GH) of the resulting model, at date ! DATE, are computed by linearly interpolating between the ! coefficients of the earlier model (GH1), at date DTE1, ! and those of the later model (GH2), at date DTE2. If one ! model is smaller than the other, the interpolation is ! performed with the missing coefficients assumed to be 0. ! !### Author ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 subroutine intershc(date, dte1, nmax1, gh1, dte2, nmax2, gh2, nmax, gh) real(wp), intent(in) :: date !! Date of resulting model (in decimal year) real(wp), intent(in) :: dte1 !! Date of earlier model integer, intent(in) :: nmax1 !! Maximum degree and order of earlier model real(wp), intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of earlier model real(wp), intent(in) :: dte2 !! Date of later model integer, intent(in) :: nmax2 !! Maximum degree and order of later model real(wp), intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of later model real(wp), intent(out) :: gh(*) !! Coefficients of resulting model integer, intent(out) :: nmax !! Maximum degree and order of resulting model real(wp) :: factor integer :: i, k, l factor = (date - dte1) / (dte2 - dte1) if (nmax1 == nmax2) then k = nmax1 * (nmax1 + 2) nmax = nmax1 elseif (nmax1 > nmax2) then k = nmax2 * (nmax2 + 2) l = nmax1 * (nmax1 + 2) do i = k + 1, l gh(i) = gh1(i) + factor * (-gh1(i)) end do nmax = nmax1 else k = nmax1 * (nmax1 + 2) l = nmax2 * (nmax2 + 2) do i = k + 1, l gh(i) = factor * gh2(i) end do nmax = nmax2 end if do i = 1, k gh(i) = gh1(i) + factor * (gh2(i) - gh1(i)) end do end subroutine intershc !***************************************************************************************** !> ! Extrapolates linearly a spherical harmonic model with a ! rate-of-change model. ! ! The coefficients (GH) of the resulting model, at date ! DATE, are computed by linearly extrapolating the coef- ! ficients of the base model (GH1), at date DTE1, using ! those of the rate-of-change model (GH2), at date DTE2. If ! one model is smaller than the other, the extrapolation is ! performed with the missing coefficients assumed to be 0. ! !### Author ! * Version 1.01, A. Zunde ! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 subroutine extrashc(date, dte1, nmax1, gh1, nmax2, gh2, nmax, gh) real(wp), intent(in) :: date !! Date of resulting model (in decimal year) real(wp), intent(in) :: dte1 !! Date of base model integer, intent(in) :: nmax1 !! Maximum degree and order of base model real(wp), intent(in) :: gh1(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of base model integer, intent(in) :: nmax2 !! Maximum degree and order of rate-of-change model real(wp), intent(in) :: gh2(*) !! Schmidt quasi-normal internal spherical harmonic coefficients of rate-of-change model real(wp), intent(out) :: gh(*) !! Coefficients of resulting model integer, intent(out) :: nmax !! Maximum degree and order of resulting model real(wp) :: factor integer :: i, k, l factor = (date - dte1) if (nmax1 == nmax2) then k = nmax1 * (nmax1 + 2) nmax = nmax1 elseif (nmax1 > nmax2) then k = nmax2 * (nmax2 + 2) l = nmax1 * (nmax1 + 2) do i = k + 1, l gh(i) = gh1(i) end do nmax = nmax1 else k = nmax1 * (nmax1 + 2) l = nmax2 * (nmax2 + 2) do i = k + 1, l gh(i) = factor * gh2(i) end do nmax = nmax2 end if do i = 1, k gh(i) = gh1(i) + factor * gh2(i) end do end subroutine extrashc !***************************************************************************************** !> ! geodetic to scaled cartesian coordinates pure function geo_to_cart(glat, glon, alt) result(x) real(wp), intent(in) :: glat !! geodetic latitude in degrees (north) real(wp), intent(in) :: glon !! geodetic longitude in degrees (east) real(wp), intent(in) :: alt !! altitude in km above sea level real(wp), dimension(3) :: x !! cartesian coordinates in earth radii (6371.2 km) !! !! * x-axis pointing to equator at 0 longitude !! * y-axis pointing to equator at 90 long. !! * z-axis pointing to north pole real(wp) :: rlat !! latitude in radians real(wp) :: rlon !! longitude in radians real(wp) :: d, rho ! deg to radians: rlat = glat * umr rlon = glon * umr ! JW : it's weird that ct is sin, and st is cos...it was like that in the original code associate (ct => sin(rlat), st => cos(rlat), cp => cos(rlon), sp => sin(rlon)) d = sqrt(aquad - (aquad - bquad) * ct * ct) rho = (alt + aquad / d) * st / era x = [rho * cp, rho * sp, (alt + bquad / d) * ct / era] end associate end function geo_to_cart end module SHELLIG_module