LCOV - code coverage report
Current view: | top level - chemistry/utils - apex.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 840 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 28 | 0.0 %

              Line data    Source code
           1             : module apex
           2             : !
           3             : ! April, 2013: B. Foster (NCAR/HAO)
           4             : !
           5             : ! This is a refactored version of the legacy apex code, originally written
           6             : ! by Art Richmond and Roy Barnes, and others in the 1995-2000 timeframe.
           7             : ! This new version is written in free-format fortran90. Subroutines and
           8             : ! module data may be use-associated from this module.
           9             : !
          10             : ! Original reference for the legacy code:
          11             : !          Richmond, A. D., Ionospheric Electrodynamics Using Magnetic Apex
          12             : !          Coordinates, J. Geomag. Geoelectr., 47, 191-212, 1995.
          13             : !
          14             : ! This code should produce near-identical results as the legacy code, altho
          15             : ! the refactored version does not provide all subroutines and options available
          16             : ! in the old code, notably the ability to write and read-back an external file.
          17             : !
          18             : ! A typical calling sequence for a code calling this module is as follows:
          19             : !
          20             : ! subroutine ggrid (legacy SUBROUTINE GGRID):
          21             : !   Make a global lat,lon,alt grid for use in later calls (optional)
          22             : !
          23             : ! subroutine apex_mka (legacy SUBROUTINE APXMKA):
          24             : !   Make magnetic arrays x,y,z,v for use in later routines
          25             : !   (geographic lat,lon grid and altitudes are input)
          26             : !   (This must be called before apex_mall and apex_q2g)
          27             : !
          28             : ! subroutine apex_mall (legacy ENTRY APXMALL):
          29             : !   Calculate modified Apex coordinates and other magnetic field parameters
          30             : !   (usually called from lat,lon,alt nested loop)
          31             : !
          32             : ! subroutine apex_q2g (legacy ENTRY APXQ2G):
          33             : !   Convert from quasi-dipole to geodetic coordinates
          34             : !   (usually called from lat,lon,alt nested loop)
          35             : !
          36             :   use shr_kind_mod,  only : r8 => shr_kind_r8
          37             :   use cam_logfile,   only : iulog
          38             :   use cam_abortutils,only : endrun
          39             :   use spmd_utils,    only : masterproc
          40             :
          41             :   implicit none
          42             :
          43             :   private
          44             :   public :: apex_set_igrf
          45             :   public :: apex_mka
          46             :   public :: apex_mall
          47             :   public :: apex_dypol
          48             :   public :: apex_subsol
          49             :   public :: apex_magloctm
          50             :   public :: apex_beg_yr
          51             :   public :: apex_end_yr
          52             :   public :: apex_q2g
          53             :
          54             :   real(r8),parameter :: re = 6371.2_r8, eps = 1.e-5_r8
          55             :
          56             :   real(r8),allocatable,save :: &
          57             :     xarray(:,:,:), & ! cos(quasi-dipole latitude)*cos(apex longitude)
          58             :     yarray(:,:,:), & ! cos(quasi-dipole latitude)*sin(apex longitude)
          59             :     zarray(:,:,:), & ! sin(quasi-dipole latitude)
          60             :     varray(:,:,:)    ! (VMP/VP)*((RE+ALT)/RE)**2
          61             : !
          62             : ! This grid (geolat,geolon,geoalt is equivalent to gdlat,gdlon,gdalt,
          63             : ! as passed to apex_mka.
          64             : !
          65             :   integer :: nglat,nglon,ngalt
          66             :   real(r8),allocatable,save :: geolat(:), geolon(:), geoalt(:)
          67             :
          68             :   integer,parameter :: nmax=13
          69             :   integer,parameter :: ncoef = nmax*nmax + 2*nmax + 1 ! 196
          70             :   real(r8),dimension(ncoef) :: &
          71             :     gb, & ! Coefficients for magnetic field calculation
          72             :     gv    ! Coefficients for magnetic potential calculation
          73             : !
          74             :   real(r8) :: &
          75             :     rtd,  & ! radians to degrees
          76             :     dtr,  & ! degrees to radians
          77             :     pola    ! pole angle (deg); when geographic lat is poleward of pola,
          78             :             ! x,y,z,v arrays are forced to be constant (pola=89.995)
          79             :
          80             :   real(r8),parameter ::       & ! Formerly common /APXCON/
          81             :     req  = 6378.160_r8,       & ! Equatorial earth radius
          82             :     precise = 7.6e-11_r8,     & ! Precision factor
          83             :     glatlim = 89.9_r8,        & ! Limit above which gradients are recalculated
          84             :     xmiss = -32767._r8
          85             : !
          86             : ! colat,elon,vp,ctp,stp were in two commons in legacy code:
          87             : ! /APXDIPL/ and /DIPOLE/. Need to check if these need to be separated.
          88             : !
          89             :   real(r8) ::  & ! Formerly /APXDIPL/ and /DIPOLE/
          90             :     colat, & ! Geocentric colatitude of geomagnetic dipole north pole (deg)
          91             :     elon,  & ! East longitude of geomagnetic dipole north pole (deg)
          92             :     vp,    & ! Magnitude, in T.m, of dipole component of magnetic
          93             :              ! potential at geomagnetic pole and geocentric radius re
          94             :     ctp,stp
          95             : !
          96             :   real(r8) ::  & ! Formerly /FLDCOMD/
          97             :     bx,    & ! X comp. of field vector at the current tracing point (Gauss)
          98             :     by,    & ! Y comp. of field vector at the current tracing point (Gauss)
          99             :     bz,    & ! Z comp. of field vector at the current tracing point (Gauss)
         100             :     bb       ! Magnitude of field vector at the current tracing point (Gauss)
         101             :
         102             :   real(r8) ::      & ! Formerly /APXIN/
         103             :     yapx(3,3)    ! Matrix of cartesian coordinates (loaded columnwise)
         104             : !
         105             : ! /ITRA/ was only in subs linapx and itrace, so can probably be removed from module data
         106             : !
         107             :   integer ::   & ! Formerly /ITRA/
         108             :     nstp         ! Step count. Incremented in sub linapx.
         109             :   real(r8)    ::   &
         110             :     y(3),      & ! Array containing current tracing point cartesian coordinates.
         111             :     yp(3),     & ! Array containing previous tracing point cartesian coordinates.
         112             :     sgn,       & ! Determines direction of trace. Set in subprogram linapx
         113             :     ds           ! Step size (Km) Computed in subprogram linapx.
         114             :
         115             :   real(r8) ::         & ! limits beyond which east-west gradients are computed
         116             :     glatmn,glatmx   ! differently to avoid potential underflow (apex_mka)
         117             :
         118             :   ! IGRF coefficients
         119             :   real(r8), allocatable :: g1(:,:), g2(:,:)
         120             :   integer :: n1, n2, ncn1, ncn2, year1, year2
         121             :   integer, protected :: apex_beg_yr
         122             :   integer, protected :: apex_end_yr
         123             :
         124             :   logical :: igrf_set = .false.
         125             :   logical :: first_warning = .false.
         126             :
         127             : contains
         128             : !-----------------------------------------------------------------------
         129             : subroutine ggrid(nvert,glatmin,glatmax,glonmin,glonmax,altmin,altmax, &
         130             :                  gplat,gplon,gpalt,mxlat,mxlon,mxalt,nlat,nlon,nalt)
         131             : !
         132             : ! Given desired range of geographic latitude, longitude and altitude,
         133             : ! choose an appropriate grid that can be used in subsequent calls to
         134             : ! subs apex_mka, apex_mall, apex_q2g.
         135             : !
         136             : ! Input args:
         137             :   integer,intent(in) :: nvert,mxlat,mxlon,mxalt
         138             :   real(r8),intent(in) :: glatmin,glatmax,glonmin,glonmax,altmin,altmax
         139             : !
         140             : ! Output args:
         141             :   integer,intent(out) :: nlat,nlon,nalt
         142             :   real(r8),intent(out) :: gplat(mxlat),gplon(mxlon),gpalt(mxalt)
         143             : !
         144             : ! Local:
         145             :   real(r8) :: dlon,dlat,diht,dnv,glonmaxx,x
         146             :   integer :: nlatmin,nlatmax,nlonmin,nlonmax,naltmin,naltmax
         147             :   integer :: i,j,k,kk
         148             :   character(len=128) :: errmsg
         149             : !
         150             : ! Check inputs:
         151             :   if (glatmin > glatmax) then
         152             :     write(errmsg,"('>>> ggrid: glatmin=',f9.2,' must be <= glatmax=',f9.2)") glatmin,glatmax
         153             :     write(iulog,*) errmsg
         154             :     call endrun( trim(errmsg) )
         155             :   endif
         156             :   if (glonmin > glonmax) then
         157             :     write(errmsg,"('>>> ggrid: glonmin=',f9.2,' must be <= glonmax=',f9.2)") glonmin,glonmax
         158             :     write(iulog,*) errmsg
         159             :     call endrun( trim(errmsg) )
         160             :   endif
         161             :   if (altmin > altmax) then
         162             :     write(errmsg,"('>>> ggrid: altmin=',f9.2,' must be <= altmax=',f9.2)") altmin,altmax
         163             :     write(iulog,*) errmsg
         164             :     call endrun( trim(errmsg) )
         165             :   endif
         166             : !
         167             : ! Init outputs:
         168             :   nlat = 0 ; nlon = 0 ; nalt = 0
         169             :   gplat = 0._r8 ; gplon = 0._r8 ; gpalt = 0._r8
         170             : !
         171             :   dnv = dble(nvert)
         172             :   dlon = 360._r8 / (5._r8*dnv)
         173             :   dlat = 180._r8 / (3._r8*dnv)
         174             :   diht = 1._r8   / dnv
         175             :
         176             :   nlatmin = max(int((glatmin+90._r8)/dlat),0)
         177             :   nlatmax = min(int((glatmax+90._r8)/dlat+1._r8),3*nvert)
         178             :   nlonmin = max(int((glonmin+180._r8)/dlon),0)
         179             :
         180             :   glonmaxx = min(glonmax,glonmin+360._r8)
         181             :   nlonmax = min(int((glonmaxx+180._r8)/dlon+1._r8),10*nvert)
         182             :
         183             :   x = re/(re+altmax)/diht-eps
         184             :   naltmin = max(x,1._r8)
         185             :   naltmin = min(naltmin,nvert-1)
         186             :   x = re/(re+altmin)/diht+eps
         187             :   i = x + 1._r8
         188             :   naltmax = min(i,nvert)
         189             :
         190             :   nlat = nlatmax - nlatmin + 1
         191             :   nlon = nlonmax - nlonmin + 1
         192             :   nlon = min(nlon,5*nvert+1)
         193             :   nalt = naltmax - naltmin + 1
         194             :
         195             :   do j=1,nlat
         196             :     gplat(j) = dlat*dble(nlatmin+j-1) - 90._r8
         197             :   enddo
         198             :   do i=1,nlon
         199             :     gplon(i) = dlon*dble(nlonmin+i-1) - 180._r8
         200             :   enddo
         201             :   do k=1,nalt
         202             :     kk = naltmax - k +1
         203             :     gpalt(k) = re*(dble(nvert-kk) - eps) / (dble(kk)+eps)
         204             :   enddo
         205             :   if (gplon(nlon-1) >= glonmax) nlon = nlon-1
         206             :   gpalt(1) = max(gpalt(1),0._r8)
         207             :
         208             :   if (masterproc) then
         209             :      write(iulog,"('ggrid: nlat=',i4,' gplat=',/,(6f9.2))") nlat,gplat
         210             :      write(iulog,"('ggrid: nlon=',i4,' gplon=',/,(6f9.2))") nlon,gplon
         211             :      write(iulog,"('ggrid: nalt=',i4,' gpalt=',/,(6f9.2))") nalt,gpalt
         212             :   endif
         213             :
         214             : end subroutine ggrid
         215             :
         216             :
         217             : !-----------------------------------------------------------------------
         218           0 : subroutine apex_set_igrf(coefs_file)
         219             :   use ioFileMod,     only : getfil
         220             :   use cam_pio_utils, only : cam_pio_openfile
         221             :   use pio,           only : file_desc_t, pio_get_var, pio_closefile, pio_nowrite, pio_inq_varid, pio_inq_dimid, pio_inq_dimlen
         222             :
         223             :   character(len=*), intent(in) :: coefs_file
         224             :
         225             :   integer :: ierr
         226             :   integer :: dim_id, var_id
         227             :   type(file_desc_t) :: ncid
         228             :   character(len=256) :: locfn
         229             :
         230           0 :   if (igrf_set) return
         231             :
         232             :   !----------------------------------------------------------------------
         233             :   !     ... open the netcdf file
         234             :   !----------------------------------------------------------------------
         235           0 :   call getfil(coefs_file, locfn, 0)
         236           0 :   call cam_pio_openfile( ncid, trim(locfn), PIO_NOWRITE)
         237             :
         238             :   !----------------------------------------------------------------------
         239             :   !     ... read the snoe dimensions
         240             :   !----------------------------------------------------------------------
         241           0 :   ierr = pio_inq_dimid( ncid, 'n1', dim_id )
         242           0 :   ierr = pio_inq_dimlen( ncid, dim_id, n1 )
         243           0 :   ierr = pio_inq_dimid( ncid, 'ncn1', dim_id )
         244           0 :   ierr = pio_inq_dimlen( ncid, dim_id, ncn1 )
         245             :
         246           0 :   ierr = pio_inq_dimid( ncid, 'n2', dim_id )
         247           0 :   ierr = pio_inq_dimlen( ncid, dim_id, n2 )
         248           0 :   ierr = pio_inq_dimid( ncid, 'ncn2', dim_id )
         249           0 :   ierr = pio_inq_dimlen( ncid, dim_id, ncn2 )
         250             :
         251           0 :   allocate( g1(n1,ncn1), g2(n2,ncn2) )
         252             :
         253           0 :   ierr = pio_inq_varid( ncid, 'g1', var_id )
         254           0 :   ierr = pio_get_var( ncid, var_id, g1 )
         255             :
         256           0 :   ierr = pio_inq_varid( ncid, 'g2', var_id )
         257           0 :   ierr = pio_get_var( ncid, var_id, g2 )
         258             :
         259           0 :   ierr = pio_inq_varid( ncid, 'era1_year', var_id )
         260           0 :   ierr = pio_get_var( ncid, var_id, year1 )
         261             :
         262           0 :   ierr = pio_inq_varid( ncid, 'era2_year', var_id )
         263           0 :   ierr = pio_get_var( ncid, var_id, year2 )
         264             :
         265           0 :   call pio_closefile(ncid)
         266             :
         267           0 :   apex_beg_yr = year1
         268           0 :   apex_end_yr = year2+5*(ncn2-1)
         269             :
         270           0 :   igrf_set = .true.
         271             :
         272           0 : end subroutine apex_set_igrf
         273             :
         274             : !-----------------------------------------------------------------------
         275           0 : subroutine apex_mka(date,gplat,gplon,gpalt,nlat,nlon,nalt,ier)
         276             : !
         277             : ! Given a 3d lat,lon,altitude grid, calculate x,y,z,v arrays in module
         278             : ! data above. These arrays are used later for calculating quantities
         279             : ! involving gradients of Apex coordinates, such as base vectors in the
         280             : ! Modified-Apex and Quasi-Dipole systems.
         281             : !
         282             : ! This defines module 3d data xarray,yarray,zarray,varray
         283             : !
         284             : ! Input args:
         285             :   real(r8),intent(in) :: date              ! year and fraction
         286             :   integer, intent(in) :: nlat,nlon,nalt ! dimensions of 3d grid
         287             :   real(r8),intent(inout) :: gplat(nlat),gplon(nlon),gpalt(nalt)
         288             : !
         289             : ! Output args:
         290             :   integer,intent(out) :: ier
         291             : !
         292             : ! Local:
         293             :   integer :: i,j,k,kpol,istat
         294             :   real(r8) :: reqore,rqorm1,cp,ct,st,sp,stmcpm,stmspm,ctm
         295             :   real(r8) :: aht,alat,phia,bmag,xmag,ymag,zdown,vmp ! apex_sub output
         296             :   real(r8) :: vnor,rp,reqam1,slp,clp,phiar
         297             :
         298           0 :   ier = 0
         299             : !
         300             : ! Some parts of the legacy apex code use constants to set dtr,rtd,
         301             : ! other parts use rtd=45./atan(1.), dtr=1./rtd. Differences are
         302             : ! on the order of 1.e-18 to 1.e-14. Here, the atan method is used.
         303             : !
         304             : !  rtd  = 5.72957795130823E1
         305             : !  dtr  = 1.745329251994330E-2
         306             : !
         307           0 :    rtd  = 45._r8/atan(1._r8)
         308           0 :    dtr  = 1._r8/rtd
         309             : !
         310             : ! pola:
         311             : !   Pole angle (deg); when the geographic latitude is poleward of POLA,
         312             : !   X,Y,Z,V are forced to be constant for all longitudes at each altitude.
         313             : !   This makes POLA = 89.995
         314             : !
         315           0 :   pola = 90._r8-sqrt(precise)*rtd    ! Pole angle (deg)
         316             :
         317             : !
         318             : ! Allocate 3d x,y,z,v arrays:
         319             : ! These are not deallocated by this module. They can be deallocated
         320             : ! by the calling program following the last call to the apex subs.
         321             : !
         322           0 :   if (.not.allocated(xarray)) then
         323           0 :     allocate(xarray(nlat,nlon,nalt),stat=istat)
         324           0 :     if (istat /= 0) call endrun( 'allocate xarray' )
         325           0 :     xarray = 0._r8
         326             :   endif
         327           0 :   if (.not.allocated(yarray)) then
         328           0 :     allocate(yarray(nlat,nlon,nalt),stat=istat)
         329           0 :     if (istat /= 0) call endrun( 'allocate yarray' )
         330           0 :     yarray = 0._r8
         331             :   endif
         332           0 :   if (.not.allocated(zarray)) then
         333           0 :     allocate(zarray(nlat,nlon,nalt),stat=istat)
         334           0 :     if (istat /= 0) call endrun( 'allocate zarray' )
         335           0 :     zarray = 0._r8
         336             :   endif
         337           0 :   if (.not.allocated(varray)) then
         338           0 :     allocate(varray(nlat,nlon,nalt),stat=istat)
         339           0 :     if (istat /= 0) call endrun( 'allocate varray' )
         340           0 :     varray = 0._r8
         341             :   endif
         342             : !
         343             : ! Set geographic grid in module data for later reference:
         344             : ! (these also are not deallocated by this module)
         345             : !
         346           0 :   nglon=nlon ; nglat=nlat ; ngalt=nalt
         347           0 :   if (.not.allocated(geolat)) allocate(geolat(nglat),stat=istat)
         348           0 :   if (.not.allocated(geolon)) allocate(geolon(nglon),stat=istat)
         349           0 :   if (.not.allocated(geoalt)) allocate(geoalt(ngalt),stat=istat)
         350           0 :   geolat(:) = gplat(:)
         351           0 :   geolon(:) = gplon(:)
         352           0 :   geoalt(:) = gpalt(:)
         353             : !
         354             : ! Set coefficients gb,gv (module data) for requested year:
         355             : !
         356           0 :   call cofrm(date)
         357             :
         358             : ! write(iulog,"('apex_mka after cofrm: ncoef=',i4,' gb=',/,(6f12.3))") ncoef,gb
         359             : ! write(iulog,"('apex_mka after cofrm: ncoef=',i4,' gv=',/,(6f12.3))") ncoef,gv
         360             :
         361           0 :   call apex_dypol(colat,elon,vp)
         362             :
         363           0 :   ctp = cos(colat*dtr)
         364           0 :   stp = sin(colat*dtr)
         365             :
         366           0 :   reqore = req/re
         367           0 :   rqorm1 = reqore-1._r8
         368             :
         369           0 :   do j=1,nlat
         370           0 :     ct = sin(gplat(j)*dtr)
         371           0 :     st = cos(gplat(j)*dtr)
         372           0 :     kpol = 0
         373           0 :     if (abs(gplat(j)) > pola) kpol = 1
         374           0 :     do i=1,nlon
         375           0 :       if (kpol==1.and.i > 1) then
         376           0 :         xarray(j,i,:) = xarray(j,1,:)
         377           0 :         yarray(j,i,:) = yarray(j,1,:)
         378           0 :         zarray(j,i,:) = zarray(j,1,:)
         379           0 :         varray(j,i,:) = varray(j,1,:)
         380             :         cycle
         381             :       endif
         382           0 :       cp = cos((gplon(i)-elon)*dtr)
         383           0 :       sp = sin((gplon(i)-elon)*dtr)
         384             : !
         385             : !  ctm   is pseudodipole component of z
         386             : ! -ctm   is pseudodipole component of v
         387             : !  stmcpm is pseudodipole component of x
         388             : !  stmspm is pseudodipole component of y
         389             : !
         390           0 :       ctm = ctp*ct + stp*st*cp
         391           0 :       stmcpm = st*ctp*cp - ct*stp
         392           0 :       stmspm = st*sp
         393           0 :       do k=1,nalt
         394           0 :         call apex_sub(date,gplat(j),gplon(i),gpalt(k),&
         395           0 :           aht,alat,phia,bmag,xmag,ymag,zdown,vmp)
         396             :
         397           0 :         vnor = vmp/vp
         398           0 :         rp = 1._r8 + gpalt(k)/re
         399           0 :         varray(j,i,k) = vnor*rp*rp + ctm
         400           0 :         reqam1 = req*(aht-1._r8)
         401           0 :         slp = sqrt(max(reqam1-gpalt(k),0._r8)/(reqam1+re))
         402             : !
         403             : ! Reverse sign of slp in southern magnetic hemisphere
         404             : !
         405           0 :         if (zdown.lt.0._r8) slp = -slp
         406           0 :         clp = sqrt (rp/(reqore*aht-rqorm1))
         407           0 :         phiar = phia*dtr
         408           0 :         xarray(j,i,k) = clp*cos (phiar) - stmcpm
         409           0 :         yarray(j,i,k) = clp*sin (phiar) - stmspm
         410           0 :         zarray(j,i,k) = slp - ctm
         411             :       enddo ! k=1,nalt
         412             :     enddo ! i=1,nlon
         413             :   enddo ! j=1,nlat
         414             : !
         415             : ! Establish for this grid polar latitude limits beyond which east-west
         416             : ! gradients are computed differently to avoid potential underflow
         417             : ! (glatmx,glatmn are in module data, glatlim is parameter constant)
         418             : !
         419           0 :   glatmx = max( glatlim,gplat(nlat-2))
         420           0 :   glatmn = min(-glatlim,gplat(2))
         421             :
         422           0 : end subroutine apex_mka
         423             : !-----------------------------------------------------------------------
         424           0 : subroutine apex_mall(glat,glon,alt,hr, b,bhat,bmag,si,alon,xlatm,vmp,w,&
         425             :   d,be3,sim,d1,d2,d3,e1,e2,e3,xlatqd,f,f1,f2,ier)
         426             : !
         427             : ! Compute Modified Apex coordinates, quasi-dipole coordinates,
         428             : ! base vectors and other parameters by interpolation from
         429             : ! precalculated arrays. Subroutine apex_mka must be called
         430             : ! before calling this subroutine.
         431             : !
         432             : ! Args:
         433             :   real(r8),intent(in)  :: & ! Input
         434             :     glat             ,& ! Geographic (geodetic) latitude (deg)
         435             :     glon             ,& ! Geographic (geodetic) longitude (deg)
         436             :     alt              ,& ! Altitude (km)
         437             :     hr                  ! Reference altitude (km)
         438             :
         439             :   real(r8),intent(out) :: & ! Output
         440             :     b(3)             ,& ! Magnetic field components (east, north, up), in nT
         441             :     bhat(3)          ,& ! components (east, north, up) of unit vector along
         442             :                         ! geomagnetic field direction
         443             :     bmag             ,& ! Magnitude of magnetic field (nT)
         444             :     si               ,& ! sin(i)
         445             :     alon             ,& ! Apex longitude = modified apex longitude =
         446             :                         ! quasi-dipole longitude (deg)
         447             :     xlatm            ,& ! Modified Apex latitude (deg)
         448             :     vmp              ,& ! Magnetic potential (T.m)
         449             :     w                ,& ! W of Richmond reference above, in km**2 /nT (i.e., 10**15 m**2 /T)
         450             :     d                ,& ! D of Richmond reference above
         451             :     be3              ,& ! B_e3 of reference above (= Bmag/D), in nT
         452             :     sim              ,& ! sin(I_m) described in Richmond reference above
         453             :     xlatqd           ,& ! Quasi-dipole latitude (deg)
         454             :     f                ,& ! F described in ref above for quasi-dipole coordinates
         455             :     f1(2),f2(2)         ! Components (east, north) of base vectors
         456             : !
         457             :   real(r8),dimension(3),intent(out) :: d1,d2,d3,e1,e2,e3 ! Components of base vectors
         458             :   integer,intent(out) :: ier ! error return
         459             : !
         460             : ! Local:
         461             :   real(r8) :: glonloc,cth,sth,glatx,clm,r3_2
         462             :   real(r8) :: fx,fy,fz,fv
         463             :   real(r8) :: dfxdth,dfydth,dfzdth,dfvdth, &
         464             :               dfxdln,dfydln,dfzdln,dfvdln, &
         465             :               dfxdh ,dfydh ,dfzdh ,dfvdh
         466             :   real(r8),dimension(3) :: gradx,grady,gradz,gradv, grclm,clmgrp,rgrlp
         467             :   real(r8) ::                    & ! dummies for polar calls to intrp
         468             :     fxdum,fydum,fzdum,fvdum,     &
         469             :     dmxdth,dmydth,dmzdth,dmvdth, &
         470             :     dmxdh,dmydh,dmzdh,dmvdh
         471             : !
         472             : ! Init:
         473             : !
         474           0 :   ier = 0
         475           0 :   glonloc = glon
         476             :
         477             :   call intrp (glat,glonloc,alt, geolat,geolon,geoalt,nglat,nglon,ngalt, &
         478             :              fx,fy,fz,fv,                                               &
         479             :              dfxdth,dfydth,dfzdth,dfvdth,                               &
         480             :              dfxdln,dfydln,dfzdln,dfvdln,                               &
         481           0 :              dfxdh ,dfydh ,dfzdh ,dfvdh, ier)
         482             :
         483           0 :   if (ier /= 0) then
         484             :     call setmiss(xmiss,xlatm,alon,vmp,b,bmag,be3,sim,si,f,d,w, &
         485           0 :       bhat,d1,d2,d3,e1,e2,e3,f1,f2)
         486             :     write(iulog,"('apex_mall called setmiss: glat,glon,alt=',3f12.3)") &
         487           0 :       glat,glon,alt
         488           0 :     return
         489             :   endif
         490             :
         491             :   call adpl(glat,glonloc,cth,sth,fx,fy,fz,fv, &
         492           0 :     dfxdth,dfydth,dfzdth,dfvdth,dfxdln,dfydln,dfzdln,dfvdln)
         493             :
         494             :   call gradxyzv(alt,cth,sth, &
         495             :     dfxdth,dfydth,dfzdth,dfvdth,dfxdln,dfydln,dfzdln,dfvdln, &
         496           0 :     dfxdh,dfydh,dfzdh,dfvdh,gradx,grady,gradz,gradv)
         497             : !
         498             : ! If the point is very close to either the North or South
         499             : ! geographic pole, recompute the east-west gradients after
         500             : ! stepping a small distance from the pole.
         501             : !
         502           0 :   if (glat > glatmx .or. glat < glatmn) then
         503           0 :     glatx = glatmx
         504           0 :     if (glat < 0._r8) glatx = glatmn
         505             :
         506             :     call intrp (glatx,glonloc,alt, geolat,geolon,geoalt,nglat,nglon,ngalt, &
         507             :                fxdum,fydum,fzdum,fvdum,                                    &
         508             :                dmxdth,dmydth,dmzdth,dmvdth,dfxdln,dfydln,dfzdln,           &
         509           0 :                dfvdln,dmxdh,dmydh,dmzdh,dmvdh, ier)
         510             :
         511             :     call adpl(glatx,glonloc,cth,sth,fxdum,fydum,fzdum,fvdum, &
         512           0 :       dmxdth,dmydth,dmzdth,dmvdth,dfxdln,dfydln,dfzdln,dfvdln)
         513             :
         514             :     call grapxyzv(alt,cth,sth,dfxdln,dfydln,dfzdln,dfvdln, &
         515           0 :       gradx,grady,gradz,gradv)
         516             :   endif
         517             :
         518             :   call gradlpv(hr,alt,fx,fy,fz,fv,gradx,grady,gradz,gradv, &
         519           0 :     xlatm,alon,vmp,grclm,clmgrp,xlatqd,rgrlp,b,clm,r3_2)
         520             :
         521             :   call basevec(hr,xlatm,grclm,clmgrp,rgrlp,b,clm,r3_2, &
         522           0 :                bmag,sim,si,f,d,w,bhat,d1,d2,d3,e1,e2,e3,f1,f2)
         523             :
         524           0 :   be3 = bmag/d
         525           0 :   ier = 0
         526             :
         527             : end subroutine apex_mall
         528             : !-----------------------------------------------------------------------
         529           0 : subroutine apex_q2g(qdlat,qdlon,alt,gdlat,gdlon,ier)
         530             : !
         531             : ! Convert from quasi-dipole to geodetic coordinates. This subroutine
         532             : ! (input magnetic, output geodetic) is the functional inverse of
         533             : ! subroutine apex_mall (input geodetic, output magnetic). Sub apex_mka
         534             : ! must be called before this routine.
         535             : !
         536             : ! Args:
         537             :   real(r8),intent(in) ::  & ! inputs
         538             :     qdlat,                & ! quasi-dipole latitude (deg)
         539             :     qdlon,                & ! quasi-dipole longitude (deg)
         540             :     alt                     ! altitude (km)
         541             :
         542             :   real(r8),intent(out) :: & ! outputs
         543             :     gdlat,            & ! geodetic latitude (deg)
         544             :     gdlon               ! geodetic longitude (deg)
         545             :   integer,intent(out) :: ier ! error return
         546             : !
         547             : ! Local:
         548             :   real(r8) :: x0,y0,z0,xnorm,xdif,ydif,zdif,dist2,hgrd2e,hgrd2n,hgrd2,&
         549             :     angdist,distlon,glatx,cal,sal,coslm,slm,cad,sad,slp,clm2,slm2,&
         550             :     sad2,cal2,clp2,clp,dylon
         551             :   real(r8) :: ylat,ylon ! first guess output by gm2gc, input to intrp
         552             :   integer :: iter
         553             :   integer,parameter :: niter=20
         554             :   real(r8) ::                    & ! output of sub intrp
         555             :     fx,fy,fz,fv,                 & ! interpolated values of x,y,z,v
         556             :     dfxdth,dfydth,dfzdth,dfvdth, & ! derivatives of x,y,z,v wrt colatitude
         557             :     dfxdln,dfydln,dfzdln,dfvdln, & ! derivatives of x,y,z,v wrt longitude
         558             :     dfxdh ,dfydh ,dfzdh ,dfvdh     ! derivatives of x,y,z,v wrt altitude
         559             :   real(r8) ::                    & ! dummies for polar calls to intrp
         560             :     fxdum,fydum,fzdum,fvdum,     &
         561             :     dmxdth,dmydth,dmzdth,dmvdth, &
         562             :     dmxdh,dmydh,dmzdh,dmvdh
         563             :   real(r8) :: cth,sth  ! output of adpl
         564             :   character(len=5) :: edge
         565             :
         566           0 :   ier = 0 ; gdlat = 0._r8 ; gdlon = 0._r8
         567             : !
         568             : ! Determine quasi-cartesian coordinates on a unit sphere of the
         569             : ! desired magnetic lat,lon in quasi-dipole coordinates.
         570             : !
         571           0 :   x0 = cos (qdlat*dtr) * cos (qdlon*dtr)
         572           0 :   y0 = cos (qdlat*dtr) * sin (qdlon*dtr)
         573           0 :   z0 = sin (qdlat*dtr)
         574             : !
         575             : ! Initial guess:  use centered dipole, convert to geocentric coords
         576             : !
         577           0 :   call gm2gc (qdlat,qdlon,ylat,ylon)
         578             : !
         579             : ! Iterate until (angular distance)**2 (units: radians) is within
         580             : ! precise of location (qdlat,qdlon) on a unit sphere.
         581             : ! (precise is a parameter in module data)
         582             : !
         583           0 :   do iter=1,niter
         584             : !
         585             : ! geolat,lon,alt and nglat,lon,alt are in module data (set by apex_mka)
         586             : !
         587             :     call intrp (ylat,ylon,alt, geolat,geolon,geoalt,nglat,nglon,ngalt, &
         588             :                fx,fy,fz,fv,                                            &
         589             :                dfxdth,dfydth,dfzdth,dfvdth,                            &
         590             :                dfxdln,dfydln,dfzdln,dfvdln,                            &
         591           0 :                dfxdh ,dfydh ,dfzdh ,dfvdh, ier)
         592           0 :     if (ier /= 0) then
         593           0 :       write(iulog,"('>>> apex_q2g error from intrp')")
         594           0 :       call endrun( 'qpex_q2g intrp'  )
         595             :     endif
         596             : !
         597             : !  Add-back of pseudodipole component to x,y,z,v and their derivatives.
         598             : !
         599             :     call adpl(ylat,ylon,cth,sth,fx,fy,fz,fv, &
         600           0 :       dfxdth,dfydth,dfzdth,dfvdth,dfxdln,dfydln,dfzdln,dfvdln)
         601           0 :     distlon = cos(ylat*dtr)
         602             :
         603           0 :     if (ylat > glatmx .or. ylat < glatmn) then ! glatmx,glatmn are module data
         604           0 :       glatx = glatmx
         605           0 :       if (ylat.lt.0._r8) glatx = glatmn
         606           0 :       distlon = cos (glatx*dtr)
         607             :       call intrp (glatx,ylon,alt, geolat,geolon,geoalt,nglat,nglon,ngalt, &
         608             :                  fxdum,fydum,fzdum,fvdum,                                 &
         609             :                  dmxdth,dmydth,dmzdth,dmvdth,dfxdln,dfydln,dfzdln,        &
         610           0 :                  dfvdln,dmxdh,dmydh,dmzdh,dmvdh, ier)
         611           0 :       if (ier /= 0) then
         612           0 :         write(iulog,"('>>> apex_q2g error from polar intrp')")
         613           0 :         call endrun( 'qpex_q2g intrp' )
         614             :       endif
         615             :
         616             :       call adpl(glatx,ylon,cth,sth,fxdum,fydum,fzdum,fvdum, &
         617           0 :         dmxdth,dmydth,dmzdth,dmvdth,dfxdln,dfydln,dfzdln,dfvdln)
         618             :     endif
         619             : !
         620             : ! At this point, FX,FY,FZ are approximate quasi-cartesian
         621             : ! coordinates on a unit sphere for the quasi-dipole coordinates
         622             : ! corresponding to the geodetic coordinates YLAT, YLON.
         623             : ! Normalize the vector length of (FX,FY,FZ) to unity using XNORM
         624             : ! so that the resultant vector can be directly compared with the
         625             : ! target vector (X0,Y0,Z0).
         626             : !
         627           0 :     xnorm = sqrt(fx*fx + fy*fy + fz*fz)
         628           0 :     xdif = fx/xnorm - x0
         629           0 :     ydif = fy/xnorm - y0
         630           0 :     zdif = fz/xnorm - z0
         631             : !
         632             : ! dist2 = square of distance between normalized (fx,fy,fz) and x0,y0,z0.
         633             : !
         634           0 :     dist2 = xdif*xdif + ydif*ydif + zdif*zdif
         635             :
         636           0 :     if (dist2 <= precise) then
         637           0 :       ier = 0
         638           0 :       gdlat = ylat
         639           0 :       gdlon = ylon
         640           0 :       return
         641             :     endif
         642             : !
         643             : ! hgrd2* = one-half of east or north gradient of dist2 on unit sphere.
         644             : !
         645           0 :     hgrd2e =  (xdif*dfxdln + ydif*dfydln + zdif*dfzdln)/distlon
         646           0 :     hgrd2n = -(xdif*dfxdth + ydif*dfydth + zdif*dfzdth)
         647           0 :     hgrd2  = sqrt(hgrd2e*hgrd2e + hgrd2n*hgrd2n)
         648             : !
         649             : ! angdist = magnitude of angular distance to be moved for new guess
         650             : !           of ylat, ylon.
         651             : !
         652           0 :     angdist = dist2/hgrd2
         653             : !
         654             : ! Following spherical trigonometry moves ylat,ylon to new location,
         655             : ! in direction of grad(dist2), by amount angdist.
         656             : !
         657           0 :     cal = -hgrd2n/hgrd2
         658           0 :     sal = -hgrd2e/hgrd2
         659           0 :     coslm = cos(ylat*dtr)
         660           0 :     slm = sin(ylat*dtr)
         661           0 :     cad = cos(angdist)
         662           0 :     sad = sin(angdist)
         663           0 :     slp = slm*cad + coslm*sad*cal
         664             :
         665           0 :     clm2 = coslm*coslm
         666           0 :     slm2 = slm*slm
         667           0 :     sad2 = sad*sad
         668           0 :     cal2 = cal*caL
         669           0 :     clp2 = clm2 + slm2*sad2 - 2._r8*slm*cad*coslm*sad*cal -clm2*sad2*cal2
         670           0 :     clp = sqrt (max(0._r8,clp2))
         671           0 :     ylat = atan2(slp,clp)*rtd
         672             : !
         673             : ! Restrict latitude iterations to stay within the interpolation grid
         674             : ! limits, but let intrp find any longitude exceedence.  This is only
         675             : ! an issue when the interpolation grid does not cover the entire
         676             : ! magnetic pole region.
         677             : !
         678           0 :     ylat = min(ylat,geolat(nglat))
         679           0 :     ylat = max(ylat,geolat(1))
         680           0 :     dylon = atan2 (sad*sal,cad*coslm-sad*slm*cal)*rtd
         681           0 :     ylon = ylon + dylon
         682           0 :     if (ylon > geolon(nglon)) ylon = ylon - 360._r8
         683           0 :     if (ylon < geolon(1))     ylon = ylon + 360._r8
         684             :
         685             :   enddo ! iter=1,niter
         686             :
         687           0 :   write(iulog,"('>>> apex_q2g: ',i3,' iterations only reduced the angular')") niter
         688             :   write(iulog,"('              difference to ',f10.5,' degrees, where test criterion')") &
         689           0 :     sqrt(dist2)*rtd
         690           0 :   write(iulog,"('              is ',f10.5,' degrees.')") sqrt(precise)*rtd
         691           0 :   edge = '     '
         692           0 :   if (ylat == geolat(nglat)) edge = 'north'
         693           0 :   if (ylat == geolat(1))     edge = 'south'
         694           0 :   if (edge /= '     ') then
         695           0 :     write(iulog,"('Coordinates are on the ',a,' edge of the interpolation grid ')") edge
         696           0 :     write(iulog,"('and latitude is constrained to stay within grid limits when iterating.')")
         697             :   endif
         698           0 :   ier = 1
         699             :
         700             : end subroutine apex_q2g
         701             : !-----------------------------------------------------------------------
         702           0 : subroutine gradxyzv(alt,cth,sth, &
         703             :     dfxdth,dfydth,dfzdth,dfvdth,dfxdln,dfydln,dfzdln,dfvdln, &
         704             :     dfxdh,dfydh,dfzdh,dfvdh,gradx,grady,gradz,gradv)
         705             : !
         706             : ! Calculates east,north,up components of gradients of x,y,z,v in
         707             : ! geodetic coordinates.  All gradients are in inverse km.  Assumes
         708             : ! flatness of 1/298.25 and equatorial radius (REQ) of 6378.16 km.
         709             : ! 940803 A. D. Richmond
         710             : !
         711             : ! Args:
         712             :   real(r8),intent(in) :: alt,cth,sth
         713             :   real(r8),dimension(3),intent(out) :: gradx,grady,gradz,gradv
         714             :   real(r8),intent(in) ::         &
         715             :     dfxdth,dfydth,dfzdth,dfvdth, &
         716             :     dfxdln,dfydln,dfzdln,dfvdln, &
         717             :     dfxdh,dfydh,dfzdh,dfvdh
         718             : !
         719             : ! Local:
         720             :   real(r8) :: d,d2,rho,dddthod,drhodth,dzetdth,ddisdth
         721             :
         722             : !
         723             : !          40680925. = req**2 (rounded off)
         724             : !          272340.   = req**2 * E2, where E2 = (2. - 1./298.25)/298.25
         725             : !                      is square of eccentricity of ellipsoid.
         726             : !
         727           0 :   d2 = 40680925.e0_r8 - 272340.e0_r8*cth*cth
         728           0 :   d = sqrt(d2)
         729           0 :   rho = sth*(alt + 40680925.e0_r8/d)
         730           0 :   dddthod = 272340.e0_r8*cth*sth/d2
         731           0 :   drhodth = alt*cth + (40680925.e0_r8/d)*(cth-sth*dddthod)
         732           0 :   dzetdth =-alt*sth - (40408585.e0_r8/d)*(sth+cth*dddthod)
         733           0 :   ddisdth = sqrt(drhodth*drhodth + dzetdth*dzetdth)
         734             :
         735           0 :   gradx(1) = dfxdln/rho
         736           0 :   grady(1) = dfydln/rho
         737           0 :   gradz(1) = dfzdln/rho
         738           0 :   gradv(1) = dfvdln/rho
         739             :
         740           0 :   gradx(2) = -dfxdth/ddisdth
         741           0 :   grady(2) = -dfydth/ddisdth
         742           0 :   gradz(2) = -dfzdth/ddisdth
         743           0 :   gradv(2) = -dfvdth/ddisdth
         744             :
         745           0 :   gradx(3) = dfxdh
         746           0 :   grady(3) = dfydh
         747           0 :   gradz(3) = dfzdh
         748           0 :   gradv(3) = dfvdh
         749             :
         750           0 : end subroutine gradxyzv
         751             : !-----------------------------------------------------------------------
         752           0 : subroutine grapxyzv(alt,cth,sth, &
         753             :   dfxdln,dfydln,dfzdln,dfvdln,gradx,grady,gradz,gradv)
         754             : !
         755             : ! Calculates east component of gradient near pole.
         756             : !
         757             : ! Args:
         758             :   real(r8),intent(in) :: alt,cth,sth
         759             :   real(r8),intent(in) :: dfxdln,dfydln,dfzdln,dfvdln
         760             :   real(r8),dimension(3),intent(inout) :: gradx,grady,gradz,gradv
         761             : !
         762             : ! Local:
         763             :   real(r8) :: d,d2,rho
         764             : !
         765             : !          40680925. = req**2 (rounded off)
         766             : !          272340.   = req**2 * E2, where E2 = (2. - 1./298.25)/298.25
         767             : !                      is square of eccentricity of ellipsoid.
         768             : !
         769           0 :   d2 = 40680925.e0_r8 - 272340.e0_r8*cth*cth
         770           0 :   d = sqrt(d2)
         771           0 :   rho = sth*(alt + 40680925.e0_r8/d)
         772             :
         773           0 :   gradx(1) = dfxdln/rho
         774           0 :   grady(1) = dfydln/rho
         775           0 :   gradz(1) = dfzdln/rho
         776           0 :   gradv(1) = dfvdln/rho
         777             :
         778           0 : end subroutine grapxyzv
         779             : !-----------------------------------------------------------------------
         780           0 : subroutine gradlpv(hr,alt,fx,fy,fz,fv,gradx,grady,gradz,gradv, &
         781             :   xlatm,xlonm,vmp,grclm,clmgrp,qdlat,rgrlp,b,clm,r3_2)
         782             : !
         783             : ! Uses gradients of x,y,z,v to compute geomagnetic field and
         784             : ! gradients of apex latitude, longitude.
         785             : !
         786             : ! Args:
         787             :   real(r8),intent(in) :: & ! scalar inputs
         788             :     hr,                  & ! reference altitude (km)
         789             :     alt,                 & ! altitude (km)
         790             :     fx,fy,fz,fv            ! interpolated values of x,y,z,v, plus
         791             :                            ! pseudodipole component
         792             :   real(r8),dimension(3),intent(in) :: & ! 3-component inputs
         793             :     gradx,grady,gradz,gradv ! interpolated gradients of x,y,z,v,
         794             :                             ! including pseudodipole components (east,north,up)
         795             : !
         796             : ! Local:
         797             :   integer :: i
         798             :   real(r8) :: rr,r,rn,sqrror,cpm,spm,bo,rn2,x2py2,xlp,slp,clp,grclp
         799             :
         800             :   real(r8),intent(out) :: & ! scalar outputs
         801             :     xlatm,  & !  modified apex latitude (lambda_m), degrees
         802             :     xlonm,  & !  apex longitude (phi_a), degrees
         803             :     vmp,    & !  magnetic potential, in T.m.
         804             :     qdlat,  & !  quasi-dipole latitude, degrees
         805             :     clm,    & !  cos(lambda_m)
         806             :     r3_2      !  ((re + alt)/(re + hr))**(3/2)
         807             :
         808             :   real(r8),dimension(3),intent(out) :: & ! 3-component outputs
         809             :     grclm,   & ! grad(cos(lambda_m)), in km-1
         810             :     clmgrp,  & ! cos(lambda_m)*grad(phi_a), in km-1
         811             :     rgrlp,   & ! (re + alt)*grad(lambda')
         812             :     b          ! magnetic field, in nT
         813             :
         814           0 :   xlatm=0._r8 ; xlonm=0._r8 ; vmp=0._r8 ; grclm=0._r8 ; clmgrp=0._r8
         815           0 :   rgrlp = 0._r8 ; b=0._r8 ; clm=0._r8 ; r3_2=0._r8 ; qdlat=0._r8
         816             :
         817           0 :   rr = re + hr
         818           0 :   r  = re + alt
         819           0 :   rn = r/re
         820           0 :   sqrror = sqrt(rr/r)
         821           0 :   r3_2 = 1._r8/sqrror/sqrror/sqrror
         822           0 :   xlonm = atan2(fy,fx)
         823           0 :   cpm = cos(xlonm)
         824           0 :   spm = sin(xlonm)
         825           0 :   xlonm = rtd*xlonm ! output
         826           0 :   bo = vp*1.e6_r8 ! vp is module data; 1.e6 converts T to nT and km-1 to m-1
         827           0 :   rn2 = rn*rn
         828           0 :   vmp = vp*fv/rn2   ! output
         829           0 :   b(1) = -bo*gradv(1)/rn2
         830           0 :   b(2) = -bo*gradv(2)/rn2
         831           0 :   b(3) = -bo*(gradv(3)-2._r8*fv/r)/rn2
         832             :
         833           0 :   x2py2 = fx*fx + fy*fy
         834           0 :   xlp = atan2(fz,sqrt(x2py2))
         835           0 :   slp = sin(xlp)
         836           0 :   clp = cos(xlp)
         837           0 :   qdlat = xlp*rtd   ! output
         838           0 :   clm = sqrror*clp  ! output
         839           0 :   if (clm > 1._r8) then
         840           0 :     write(iulog,"('>>> gradlpv: hr=',f12.3,' alt=',f12.3)") hr,alt
         841           0 :     write(iulog,"('    Point lies below field line that peaks at reference height.')")
         842           0 :     call endrun( 'gradlpv' )
         843             :   endif
         844           0 :   xlatm = rtd*acos(clm)
         845             : !
         846             : !  If southern magnetic hemisphere, reverse sign of xlatm
         847             : !
         848           0 :   if (slp < 0._r8) xlatm = -xlatm
         849           0 :   do i=1,3
         850           0 :     grclp = cpm*gradx(i) + spm*grady(i)
         851           0 :     rgrlp(i) = r*(clp*gradz(i) - slp*grclp)
         852           0 :     grclm(i) = sqrror*grclp
         853           0 :     clmgrp(i) = sqrror*(cpm*grady(i)-spm*gradx(i))
         854             :   enddo
         855           0 :   grclm(3) = grclm(3) - sqrror*clp/(2._r8*r)
         856             :
         857           0 : end subroutine gradlpv
         858             : !-----------------------------------------------------------------------
         859           0 : subroutine basevec(hr,xlatm,grclm,clmgrp,rgrlp,b,clm,r3_2, &
         860             :                    bmag,sim,si,f,d,w,bhat,d1,d2,d3,e1,e2,e3,f1,f2)
         861             : !
         862             : ! Computes base vectors and other parameters for apex coordinates.
         863             : ! Vector components:  east, north, up
         864             : !
         865             : ! Args:
         866             :   real(r8),intent(in) :: & ! scalar inputs
         867             :     hr,      & ! reference altitude
         868             :     xlatm,   & ! modified apex latitude (deg)
         869             :     clm,     & ! cos(lambda_m)
         870             :     r3_2       ! ((re + altitude)/(re + hr))**(3/2)
         871             :
         872             :   real(r8),dimension(3),intent(in) :: & ! 3-component inputs
         873             :     grclm,   & ! grad(cos(lambda_m)), in km-1
         874             :     clmgrp,  & ! cos(lambda_m)*grad(phi_a), in km-1
         875             :     rgrlp,   & ! (re + altitude)*grad(lambda')
         876             :     b          ! ((re + altitude)/(re + hr))**(3/2)
         877             :
         878             :   real(r8),intent(out) :: & ! scalar output
         879             :     bmag,    & ! magnitude of magnetic field, in nT
         880             :     sim,     & ! sin(I_m) of Richmond reference
         881             :     si,      & ! sin(I)
         882             :     f,       & ! F of Richmond reference
         883             :     d,       & ! D of Richmond reference
         884             :     w          ! W of Richmond reference
         885             :
         886             :   real(r8),dimension(3),intent(out) :: & ! 3-component outputs
         887             :     bhat,             & ! unit vector along geomagnetic field direction
         888             :     d1,d2,d3,e1,e2,e3   ! base vectors of Richmond reference
         889             :   real(r8),dimension(2),intent(out) :: & ! 2-component outputs
         890             :     f1,f2               ! base vectors of Richmond reference
         891             : !
         892             : ! Local:
         893             :   integer :: i
         894             :   real(r8) :: rr,simoslm,d1db,d2db
         895             :
         896           0 :   rr = re + hr
         897           0 :   simoslm = 2._r8/sqrt(4._r8 - 3._r8*clm*clm)
         898           0 :   sim = simoslm*sin(xlatm*dtr)
         899           0 :   bmag = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
         900           0 :   d1db = 0._r8
         901           0 :   d2db = 0._r8
         902           0 :   do i=1,3
         903           0 :     bhat(i) = b(i)/bmag
         904           0 :     d1(i) = rr*clmgrp(i)
         905           0 :     d1db = d1db + d1(i)*bhat(i)
         906           0 :     d2(i) = rr*simoslm*grclm(i)
         907           0 :     d2db = d2db + d2(i)*bhat(i)
         908             :   enddo
         909             : !
         910             : ! Ensure that d1,d2 are exactly perpendicular to B:
         911             : !
         912           0 :   do i=1,3
         913           0 :     d1(i) = d1(i) - d1db*bhat(i)
         914           0 :     d2(i) = d2(i) - d2db*bhat(i)
         915             :   enddo
         916           0 :   e3(1) = d1(2)*d2(3) - d1(3)*d2(2)
         917           0 :   e3(2) = d1(3)*d2(1) - d1(1)*d2(3)
         918           0 :   e3(3) = d1(1)*d2(2) - d1(2)*d2(1)
         919           0 :   d = bhat(1)*e3(1) + bhat(2)*e3(2) + bhat(3)*e3(3)
         920           0 :   do i=1,3
         921           0 :     d3(i) = bhat(i)/d
         922           0 :     e3(i) = bhat(i)*d ! Ensure that e3 lies along bhat.
         923             :   enddo
         924           0 :   e1(1) = d2(2)*d3(3) - d2(3)*d3(2)
         925           0 :   e1(2) = d2(3)*d3(1) - d2(1)*d3(3)
         926           0 :   e1(3) = d2(1)*d3(2) - d2(2)*d3(1)
         927           0 :   e2(1) = d3(2)*d1(3) - d3(3)*d1(2)
         928           0 :   e2(2) = d3(3)*d1(1) - d3(1)*d1(3)
         929           0 :   e2(3) = d3(1)*d1(2) - d3(2)*d1(1)
         930           0 :   w = rr*rr*clm*abs(sim)/(bmag*d)
         931           0 :   si = -bhat(3)
         932           0 :   f1(1) =  rgrlp(2)
         933           0 :   f1(2) = -rgrlp(1)
         934           0 :   f2(1) = -d1(2)*r3_2
         935           0 :   f2(2) =  d1(1)*r3_2
         936           0 :   f = f1(1)*f2(2) - f1(2)*f2(1)
         937             :
         938           0 : end subroutine basevec
         939             : !-----------------------------------------------------------------------
         940           0 : subroutine apex_dypol(colat,elon,vp)
         941             : !
         942             : ! Output args:
         943             :   real(r8),intent(out) :: &
         944             :     colat, & ! Geocentric colatitude of geomagnetic dipole north pole (deg)
         945             :     elon,  & ! East longitude of geomagnetic dipole north pole (deg)
         946             :     vp       ! Magnitude, in T.m, of dipole component of magnetic
         947             :              ! potential at geomagnetic pole and geocentric radius re
         948             : !
         949             : ! Local:
         950             :   real(r8) :: gpl,ctp
         951             : !
         952             : ! Compute geographic colatitude and longitude of the north pole of
         953             : ! earth centered dipole
         954             : !
         955           0 :   gpl = sqrt( gb(2  )**2+ gb(3  )**2+ gb(4  )**2)
         956           0 :   ctp = gb(2  )/gpl
         957             :
         958           0 :   colat = (acos(ctp))*rtd
         959           0 :   elon = atan2( gb(4  ), gb(3  ))*rtd
         960             : !
         961             : ! Compute magnitude of magnetic potential at pole, radius Re.
         962             : !      .2 = 2*(10**-4 T/gauss)*(1000 m/km) (2 comes through f0 in COFRM).
         963             : !
         964           0 :   vp = .2_r8*gpl*re
         965             :
         966           0 : end subroutine apex_dypol
         967             : !-----------------------------------------------------------------------
         968           0 : subroutine apex_sub(date,dlat,dlon,alt,aht,alat,alon,bmag,xmag,ymag,zmag,vmp)
         969             : !
         970             : ! Args:
         971             :   real(r8),intent(in) :: date
         972             :   real(r8),intent(inout) :: dlat,dlon,alt
         973             :   real(r8),intent(out) :: aht,alat,alon,bmag,xmag,ymag,zmag,vmp
         974             : !
         975             : ! Local:
         976             :   real(r8) :: clatp,polon,vpol,x,y,z,xre,yre,zre
         977             :   integer :: iflag
         978             :
         979           0 :   call cofrm(date)
         980           0 :   call apex_dypol(clatp,polon,vpol)
         981             : !
         982             : ! colat,ctp,stp,elon,vp are in module data.
         983             : !
         984           0 :   colat = clatp
         985           0 :   ctp   = cos(clatp*dtr)
         986           0 :   stp   = sqrt(1._r8-ctp*ctp)
         987             :
         988           0 :   elon  = polon
         989           0 :   vp    = vpol
         990             :
         991           0 :   vmp = 0._r8
         992             : !
         993             : ! Last 7 args of linapx are output:
         994             : !
         995           0 :   call linapx(dlat,dlon,alt, aht,alat,alon,xmag,ymag,zmag,bmag)
         996             :
         997           0 :   xmag = xmag*1.e5_r8
         998           0 :   ymag = ymag*1.e5_r8
         999           0 :   zmag = zmag*1.e5_r8
        1000           0 :   bmag = bmag*1.e5_r8
        1001           0 :   call gd2cart (dlat,dlon,alt,x,y,z)
        1002           0 :   iflag = 3
        1003           0 :   xre = x/re ; yre = y/re ; zre = z/re
        1004           0 :   call feldg(iflag,xre,yre,zre,bx,by,bz,vmp)
        1005             :
        1006           0 : end subroutine apex_sub
        1007             : !-----------------------------------------------------------------------
        1008           0 : subroutine linapx(gdlat,glon,alt,aht,alat,alon,xmag,ymag,zmag,fmag)
        1009             : !
        1010             : ! Input Args:
        1011             : !
        1012             :   real(r8),intent(inout) :: & ! These may be changed by convrt, depending on iflag
        1013             :     gdlat,                  & ! latitude of starting point (deg)
        1014             :     glon,                   & ! longitude of starting point (deg)
        1015             :     alt                       ! height of starting point (km)
        1016             : !
        1017             : ! Output Args:
        1018             : !
        1019             :   real(r8),intent(out) :: &
        1020             :     aht,              & ! (Apex height+req)/req, where req is equatorial earth radius
        1021             :     alat,             & ! Apex latitude (deg)
        1022             :     alon,             & ! Apex longitude (deg)
        1023             :     xmag,             & ! North component of magnetic field at starting point
        1024             :     ymag,             & ! East component of magnetic field at starting point
        1025             :     zmag,             & ! Down component of magnetic field at starting point
        1026             :     fmag                ! Magnetic field magnitude at starting point
        1027             : !
        1028             : ! Local:
        1029             : !
        1030             :   real(r8) :: gclat,r,singml,cgml2,rho,xlat,xlon,ht
        1031             :   real(r8) :: bnrth,beast,bdown,babs,y1,y2,y3
        1032             :   integer :: iflag,iapx
        1033             :   integer,parameter :: maxs = 200
        1034             : !
        1035             : ! Determine step size as a function of geomagnetic dipole
        1036             : ! coordinates of the starting point
        1037             : !
        1038           0 :   iflag = 2 ! gclat,r are returned
        1039           0 :   call convrt(iflag,gdlat,alt,gclat,r)
        1040             :
        1041           0 :   singml = ctp*sin(gclat*dtr) + stp*cos(gclat*dtr)*cos((glon-elon)*dtr)
        1042           0 :   cgml2 = max(0.25_r8,1._r8-singml*singml)
        1043           0 :   ds = .06_r8*r/cgml2 - 370._r8 ! ds is in module data
        1044             :
        1045           0 :   yapx = 0._r8 ! init (module data)
        1046             : !
        1047             : ! Convert from geodetic to earth centered cartesian coordinates:
        1048             : !
        1049           0 :   call gd2cart(gdlat,glon,alt,y(1),y(2),y(3))
        1050           0 :   nstp = 0
        1051             : !
        1052             : ! Get magnetic field components to determine the direction for tracing field line:
        1053             : !
        1054           0 :   iflag = 1
        1055           0 :   call feldg(iflag,gdlat,glon,alt,xmag,ymag,zmag,fmag)
        1056             :
        1057           0 :   sgn = sign(1._r8,-zmag)
        1058             : !
        1059             : ! Use cartesian coordinates to get magnetic field components
        1060             : ! (from which gradients steer the tracing)
        1061             : !
        1062             : 100 continue
        1063           0 :   iflag = 2 ! module data bx,by,bz,bb are returned
        1064           0 :   y1 = y(1)/re ; y2 = y(2)/re ; y3 = y(3)/re
        1065           0 :   call feldg(iflag,y1,y2,y3,bx,by,bz,bb)
        1066           0 :   nstp = nstp + 1
        1067             : !
        1068             : ! Quit if too many steps.
        1069             : !
        1070           0 :   if (nstp >= maxs) then
        1071           0 :     rho = sqrt(y(1)*y(1) + y(2)*y(2))
        1072           0 :     iflag = 3 ! xlat and ht are returned
        1073           0 :     call convrt(iflag,xlat,ht,rho,y(3))
        1074           0 :     xlon = rtd*atan2(y(2),y(1))
        1075           0 :     iflag = 1
        1076           0 :     call feldg(iflag,xlat,xlon,ht,bnrth,beast,bdown,babs)
        1077           0 :     call dipapx(xlat,xlon,ht,bnrth,beast,bdown,aht,alon)
        1078           0 :     alat = -sgn*rtd*acos(sqrt(1._r8/aht))
        1079           0 :     return
        1080             :   endif
        1081             : !
        1082             : ! Find next point using adams algorithm after 7 points
        1083             : !
        1084           0 :   call itrace(iapx)
        1085           0 :   if (iapx == 1) goto 100
        1086             : !
        1087             : ! Maximum radius just passed.  Find apex coords
        1088             : !
        1089           0 :   call fndapx(alt,zmag,aht,alat,alon)
        1090             :
        1091             : end subroutine linapx
        1092             : !-----------------------------------------------------------------------
        1093           0 : subroutine convrt(iflag,gdlat,alt,x1,x2)
        1094             : !
        1095             : ! Convert space point from geodetic to geocentric or vice versa.
        1096             : !
        1097             : ! iflag = 1: Convert from geodetic to cylindrical
        1098             : !   Input:  gdlat = Geodetic latitude (deg)
        1099             : !           alt   = Altitude above reference ellipsoid (km)
        1100             : !   Output: x1    = Distance from Earth's rotation axis (km)
        1101             : !           x2    = Distance above (north of) Earth's equatorial plane (km)
        1102             : !
        1103             : ! iflag = 2: Convert from geodetic to geocentric spherical
        1104             : !   Input:  gdlat = Geodetic latitude (deg)
        1105             : !           alt   = Altitude above reference ellipsoid (km)
        1106             : !   Output: x1    = Geocentric latitude (deg)
        1107             : !           x2    = Geocentric distance (km)
        1108             : !
        1109             : ! iflag = 3: Convert from cylindrical to geodetic
        1110             : !   Input:  x1    = Distance from Earth's rotation axis (km)
        1111             : !           x2    = Distance from Earth's equatorial plane (km)
        1112             : !   Output: gdlat = Geodetic latitude (deg)
        1113             : !           alt   = Altitude above reference ellipsoid (km)
        1114             : !
        1115             : ! iflag = 4: Convert from geocentric spherical to geodetic
        1116             : !   Input:  x1    = Geocentric latitude (deg)
        1117             : !           x2    = Geocentric distance (km)
        1118             : !   Output: gdlat = Geodetic latitude (deg)
        1119             : !           alt   = Altitude above reference ellipsoid (km)
        1120             : !
        1121             : ! Args:
        1122             :   integer,intent(in) :: iflag
        1123             :   real(r8),intent(inout) :: gdlat,alt
        1124             :   real(r8),intent(inout) :: x1,x2
        1125             : !
        1126             : ! Local:
        1127             :   real(r8) :: sinlat,coslat,d,z,rho,rkm,scl,gclat,ri,a2,a4,a6,a8,&
        1128             :     ccl,s2cl,c2cl,s4cl,c4cl,s8cl,s6cl,dltcl,sgl
        1129             :   real(r8),parameter ::                                                &
        1130             :     fltnvrs = 298.25_r8                                              , &
        1131             :     e2=(2._r8-1._r8/fltnvrs)/fltnvrs                                 , &
        1132             :     e4=e2*e2, e6=e4*e2, e8=e4*e4                                     , &
        1133             :     ome2req = (1._r8-e2)*req                                         , &
        1134             :     A21 = (512._r8*E2 + 128._r8*E4 + 60._r8*E6 + 35._r8*E8)/1024._r8 , &
        1135             :     A22 = (                                 E6 +        E8)/  32._r8 , &
        1136             :     A23 = -3._r8*(                    4._r8*E6 +  3._r8*E8)/ 256._r8 , &
        1137             :     A41 =    -(          64._r8*E4 + 48._r8*E6 + 35._r8*E8)/1024._r8 , &
        1138             :     A42 =     (           4._r8*E4 +  2._r8*E6 +        E8)/  16._r8 , &
        1139             :     A43 =                                        15._r8*E8 / 256._r8 , &
        1140             :     A44 =                                              -E8 /  16._r8 , &
        1141             :     A61 =  3._r8*(                    4._r8*E6 +  5._r8*E8)/1024._r8 , &
        1142             :     A62 = -3._r8*(                          E6 +        E8)/  32._r8 , &
        1143             :     A63 = 35._r8*(                    4._r8*E6 +  3._r8*E8)/ 768._r8 , &
        1144             :     A81 =                                        -5._r8*E8 /2048._r8 , &
        1145             :     A82 =                                        64._r8*E8 /2048._r8 , &
        1146             :     A83 =                                      -252._r8*E8 /2048._r8 , &
        1147             :     A84 =                                       320._r8*E8 /2048._r8
        1148             :
        1149           0 :   if (iflag < 3) then ! geodetic to geocentric
        1150             : !
        1151             : ! Compute rho,z
        1152           0 :     sinlat = sin(gdlat*dtr)
        1153           0 :     coslat = sqrt(1._r8-sinlat*sinlat)
        1154           0 :     d      = sqrt(1._r8-e2*sinlat*sinlat)
        1155           0 :     z      = (alt+ome2req/d)*sinlat
        1156           0 :     rho    = (alt+req/d)*coslat
        1157           0 :     x1 = rho
        1158           0 :     x2 = z
        1159           0 :     if (iflag == 1) return
        1160             : !
        1161             : ! Compute gclat,rkm
        1162           0 :     rkm   = sqrt(z*z+rho*rho)
        1163           0 :     gclat = rtd*atan2(z,rho)
        1164           0 :     x1 = gclat
        1165           0 :     x2 = rkm
        1166           0 :     return    ! iflag == 2
        1167             :   endif ! iflag < 3
        1168             : !
        1169             : ! Geocentric to geodetic:
        1170           0 :   if (iflag == 3) then
        1171           0 :     rho = x1
        1172           0 :     z   = x2
        1173           0 :     rkm = sqrt(z*z+rho*rho)
        1174           0 :     scl = z/rkm
        1175           0 :     gclat = asin(scl)*rtd
        1176           0 :   elseif (iflag == 4) then
        1177           0 :     gclat = x1
        1178           0 :     rkm = x2
        1179           0 :     scl = sin(gclat*dtr)
        1180             :   else
        1181             :     return
        1182             :   endif
        1183             : !
        1184             : ! iflag == 3 or 4:
        1185             : !
        1186           0 :   ri = req/rkm
        1187           0 :   a2 = ri*(a21+ri*(a22+ri* a23))
        1188           0 :   a4 = ri*(a41+ri*(a42+ri*(a43+ri*a44)))
        1189           0 :   a6 = ri*(a61+ri*(a62+ri* a63))
        1190           0 :   a8 = ri*(a81+ri*(a82+ri*(a83+ri*a84)))
        1191           0 :   ccl = sqrt(1._r8-scl*scl)
        1192           0 :   s2cl = 2._r8*scl*ccL
        1193           0 :   c2cl = 2._r8*ccl*ccl-1._r8
        1194           0 :   s4cl = 2._r8*s2cl*c2cl
        1195           0 :   c4cl = 2._r8*c2cl*c2cl-1._r8
        1196           0 :   s8cl = 2._r8*s4cl*c4cl
        1197           0 :   s6cl = s2cl*c4cl+c2cl*s4cl
        1198           0 :   dltcl = s2cl*a2+s4cl*a4+s6cl*a6+s8cl*a8
        1199           0 :   gdlat = dltcl*rtd+gclat
        1200           0 :   sgl = sin(gdlat*dtr)
        1201           0 :   alt = rkm*cos(dltcl)-req*sqrt(1._r8-e2*sgl*sgl)
        1202             :
        1203             : end subroutine convrt
        1204             : !-----------------------------------------------------------------------
        1205           0 : subroutine gd2cart(gdlat,glon,alt,x,y,z)
        1206             : !
        1207             : ! Arg:
        1208             :   real(r8),intent(inout) :: gdlat,alt,z
        1209             :   real(r8),intent(in) :: glon
        1210             :   real(r8),intent(out) :: x,y
        1211             : !
        1212             : ! Local:
        1213             :   real(r8) :: ang,rho
        1214             :   integer :: iflag
        1215             :
        1216           0 :   iflag = 1 ! Convert from geodetic to cylindrical (rho,z are output)
        1217           0 :   call convrt(iflag,gdlat,alt,rho,z)
        1218             :
        1219           0 :   ang = glon*dtr
        1220           0 :   x = rho*cos(ang)
        1221           0 :   y = rho*sin(ang)
        1222             :
        1223           0 : end subroutine gd2cart
        1224             : !-----------------------------------------------------------------------
        1225           0 : subroutine feldg(iflag,glat,glon,alt,bnrth,beast,bdown,babs)
        1226             : !
        1227             : ! Compute the DGRF/IGRF field components at the point glat,glon,alt.
        1228             : ! cofrm must be called to establish coefficients for correct date
        1229             : ! prior to calling FELDG.
        1230             : !
        1231             : ! iflag = 1:
        1232             : !   Inputs:
        1233             : !     glat = Latitude of point (deg)
        1234             : !     glon = Longitude of point (deg)
        1235             : !     alt  = Height of point (km)
        1236             : !   Outputs:
        1237             : !     bnrth = North component of field vector (Gauss)
        1238             : !     beast = East component of field vector (Gauss)
        1239             : !     bdown = Downward component of field vector (Gauss)
        1240             : !     babs  = Magnitude of field vector (Gauss)
        1241             : !
        1242             : ! iflag = 2:
        1243             : !   Inputs:
        1244             : !     glat = x coordinate (in units of earth radii 6371.2 km)
        1245             : !     glon = y coordinate (in units of earth radii 6371.2 km)
        1246             : !     alt  = z coordinate (in units of earth radii 6371.2 km)
        1247             : !   Outputs:
        1248             : !     bnrth = x component of field vector (Gauss)
        1249             : !     beast = y component of field vector (Gauss)
        1250             : !     bdown = z component of field vector (Gauss)
        1251             : !     babs  = Magnitude of field vector (Gauss)
        1252             : !
        1253             : ! iflag = 3:
        1254             : !   Inputs:
        1255             : !     glat = x coordinate (in units of earth radii 6371.2 km)
        1256             : !     glon = y coordinate (in units of earth radii 6371.2 km)
        1257             : !     alt  = z coordinate (in units of earth radii 6371.2 km)
        1258             : !   Outputs:
        1259             : !     bnrth = Dummy variable
        1260             : !     beast = Dummy variable
        1261             : !     babs  = Legacy code had "Dummy variable" here, but its
        1262             : !             set at the end if iflag==3.
        1263             : !
        1264             : ! Args:
        1265             :   integer,intent(in)     :: iflag
        1266             :   real(r8),intent(in)    :: glon
        1267             :   real(r8),intent(inout) :: glat
        1268             :   real(r8),intent(inout) :: alt
        1269             :   real(r8),intent(out)   :: bnrth,beast,bdown,babs
        1270             : !
        1271             : ! Local:
        1272             :   integer :: i,is,ihmax,last,imax,mk,k,ih,m,il,ihm,ilm
        1273             :   real(r8) :: rlat,ct,st,rlon,cp,sp,xxx,yyy,zzz,rq,f,x,y,z
        1274             :   real(r8) :: xi(3),h(ncoef),g(ncoef)
        1275             :   real(r8) :: s,t,bxxx,byyy,bzzz,brho
        1276             :
        1277           0 :   if (iflag == 1) then
        1278           0 :     is   = 1
        1279           0 :     rlat = glat*dtr
        1280           0 :     ct   = sin(rlat)
        1281           0 :     st   = cos(rlat)
        1282           0 :     rlon = glon*dtr
        1283           0 :     cp   = cos(rlon)
        1284           0 :     sp   = sin(rlon)
        1285           0 :     call gd2cart(glat,glon,alt,xxx,yyy,zzz)
        1286           0 :     xxx = xxx/re
        1287           0 :     yyy = yyy/re
        1288           0 :     zzz = zzz/re
        1289             :   else
        1290           0 :     is  = 2
        1291           0 :     xxx = glat
        1292           0 :     yyy = glon
        1293           0 :     zzz = alt
        1294             :   endif
        1295           0 :   rq    = 1._r8/(xxx**2+yyy**2+zzz**2)
        1296           0 :   xi(1) = xxx*rq
        1297           0 :   xi(2) = yyy*rq
        1298           0 :   xi(3) = zzz*rq
        1299           0 :   ihmax = nmax*nmax+1
        1300           0 :   last  = ihmax+nmax+nmax
        1301           0 :   imax  = nmax+nmax-1
        1302             : !
        1303             : ! Legacy code checks here to see if iflag or last call to cofrm have changed.
        1304             : ! For now, just do it anyway:
        1305             : !
        1306           0 :   if (iflag /= 3) then
        1307           0 :     do i=1,last
        1308           0 :       g(i) = gb(i) ! gb is module data from cofrm
        1309             :     enddo
        1310             :   else
        1311           0 :     do i=1,last
        1312           0 :       g(i) = gv(i) ! gv is module data from cofrm
        1313             :     enddo
        1314             :   endif
        1315             :
        1316           0 :   do i=ihmax,last
        1317           0 :     h(i) = g(i)
        1318             :   enddo
        1319             :
        1320             :   mk = 3
        1321             :   if (imax == 1) mk = 1
        1322             :
        1323           0 :   do k=1,mk,2
        1324           0 :     i = imax
        1325           0 :     ih = ihmax
        1326             :
        1327             : 100 continue
        1328           0 :     il = ih-i
        1329           0 :     f = 2._r8/dble(i-k+2)
        1330           0 :     x = xi(1)*f
        1331           0 :     y = xi(2)*f
        1332           0 :     z = xi(3)*(f+f)
        1333             :
        1334           0 :     i = i-2
        1335           0 :     if (i < 1) then
        1336           0 :       h(il) = g(il) + z*h(ih) + 2._r8*(x*h(ih+1)+y*h(ih+2))
        1337           0 :     elseif (i == 1) then
        1338           0 :       h(il+2) = g(il+2) + z*h(ih+2) + x*h(ih+4) - y*(h(ih+3)+h(ih))
        1339           0 :       h(il+1) = g(il+1) + z*h(ih+1) + y*h(ih+4) + x*(h(ih+3)-h(ih))
        1340           0 :       h(il)   = g(il)   + z*h(ih)   + 2._r8*(x*h(ih+1)+y*h(ih+2))
        1341             :     else
        1342           0 :       do m=3,i,2
        1343           0 :         ihm = ih+m
        1344           0 :         ilm = il+m
        1345           0 :         h(ilm+1) = g(ilm+1)+ z*h(ihm+1) + x*(h(ihm+3)-h(ihm-1))- &
        1346           0 :                    y*(h(ihm+2)+h(ihm-2))
        1347           0 :         h(ilm)   = g(ilm)  + z*h(ihm)   + x*(h(ihm+2)-h(ihm-2))+ &
        1348           0 :                    y*(h(ihm+3)+h(ihm-1))
        1349             :       enddo
        1350           0 :       h(il+2) = g(il+2) + z*h(ih+2) + x*h(ih+4) - y*(h(ih+3)+h(ih))
        1351           0 :       h(il+1) = g(il+1) + z*h(ih+1) + y*h(ih+4) + x*(h(ih+3)-h(ih))
        1352           0 :       h(il)   = g(il)   + z*h(ih)   + 2._r8*(x*h(ih+1)+y*h(ih+2))
        1353             :     endif
        1354             :
        1355           0 :     ih = il
        1356           0 :     if (i >= k) goto 100
        1357             :   enddo ! k=1,mk,2
        1358             :
        1359           0 :   s = .5_r8*h(1)+2._r8*(h(2)*xi(3)+h(3)*xi(1)+h(4)*xi(2))
        1360           0 :   t = (rq+rq)*sqrt(rq)
        1361           0 :   bxxx = t*(h(3)-s*xxx)
        1362           0 :   byyy = t*(h(4)-s*yyy)
        1363           0 :   bzzz = t*(h(2)-s*zzz)
        1364           0 :   babs = sqrt(bxxx**2+byyy**2+bzzz**2)
        1365           0 :   if (is .eq. 1) then            ! (convert back to geodetic)
        1366           0 :     beast = byyy*cp-bxxx*sp
        1367           0 :     brho  = byyy*sp+bxxx*cp
        1368           0 :     bnrth = bzzz*st-brho*ct
        1369           0 :     bdown = -bzzz*ct-brho*st
        1370           0 :   elseif (is .eq. 2) then        ! (leave in earth centered cartesian)
        1371           0 :     bnrth = bxxx
        1372           0 :     beast = byyy
        1373           0 :     bdown = bzzz
        1374             :   endif
        1375             : !
        1376             : ! Magnetic potential computation makes use of the fact that the
        1377             : ! calculation of V is identical to that for r*Br, if coefficients
        1378             : ! in the latter calculation have been divided by (n+1) (coefficients
        1379             : ! GV).  Factor .1 converts km to m and gauss to tesla.
        1380             : !
        1381           0 :   if (iflag == 3) babs = (bxxx*xxx + byyy*yyy + bzzz*zzz)*re*.1_r8
        1382             :
        1383           0 : end subroutine feldg
        1384             : !-----------------------------------------------------------------------
        1385           0 : subroutine dipapx(gdlat,gdlon,alt,bnorth,beast,bdown,a,alon)
        1386             : !
        1387             : ! Compute a, alon from local magnetic field using dipole and spherical approx.
        1388             : ! Reference from legacy code: 940501 A. D. Richmond
        1389             : !
        1390             : ! Input:
        1391             : !   gdlat  = geodetic latitude, degrees
        1392             : !   gdlon  = geodetic longitude, degrees
        1393             : !   alt    = altitude, km
        1394             : !   bnorth = geodetic northward magnetic field component (any units)
        1395             : !   beast  = eastward magnetic field component
        1396             : !   bdown  = geodetic downward magnetic field component
        1397             : ! Output:
        1398             : !   a      = apex radius, 1 + h_A/R_eq
        1399             : !   alon   = apex longitude, degrees
        1400             : !
        1401             : ! Algorithm:
        1402             : !   Use spherical coordinates.
        1403             : !   Let GP be geographic pole.
        1404             : !   Let GM be geomagnetic pole (colatitude COLAT, east longitude ELON).
        1405             : !   Let G be point at GDLAT,GDLON.
        1406             : !   Let E be point on sphere below apex of dipolar field line passing through G.
        1407             : !   Let TD be dipole colatitude of point G, found by applying dipole formula
        1408             : !     for dip angle to actual dip angle.
        1409             : !   Let B be Pi plus local declination angle.  B is in the direction
        1410             : !     from G to E.
        1411             : !   Let TG be colatitude of G.
        1412             : !   Let ANG be longitude angle from GM to G.
        1413             : !   Let TE be colatitude of E.
        1414             : !   Let TP be colatitude of GM.
        1415             : !   Let A be longitude angle from G to E.
        1416             : !   Let APANG = A + ANG
        1417             : !   Let PA be geomagnetic longitude, i.e., Pi minus angle measured
        1418             : !     counterclockwise from arc GM-E to arc GM-GP.
        1419             : !   Let TF be arc length between GM and E.
        1420             : !   Then, using notation C=cos, S=sin, COT=cot, spherical-trigonometry formulas
        1421             : !     for the functions of the angles are as shown below.  Note: STFCPA,
        1422             : !     STFSPA are sin(TF) times cos(PA), sin(PA), respectively.
        1423             : !
        1424             :   real(r8),intent(in)  :: gdlat,gdlon,alt,bnorth,beast,bdown
        1425             :   real(r8),intent(out) :: a,alon
        1426             : !
        1427             : ! Local:
        1428             :   real(r8) :: bhor,std,ctd,sb,cb,ctg,stg,ang,sang,cang,cte,ste,sa,ca, &
        1429             :               cottd,capang,sapang,stfcpa,stfspa,ha,r
        1430             :
        1431           0 :   bhor = sqrt(bnorth*bnorth + beast*beast)
        1432           0 :   if (bhor == 0._r8) then
        1433           0 :     alon = 0._r8
        1434           0 :     a = 1.e34_r8
        1435           0 :     return
        1436             :   endif
        1437             :
        1438           0 :   cottd = bdown*.5_r8/bhor
        1439           0 :   std = 1._r8/sqrt(1._r8+cottd*cottd)
        1440           0 :   ctd = cottd*std
        1441           0 :   sb = -beast/bhor
        1442           0 :   cb = -bnorth/bhor
        1443           0 :   ctg = sin(gdlat*dtr)
        1444           0 :   stg = cos(gdlat*dtr)
        1445           0 :   ang = (gdlon-elon)*dtr
        1446           0 :   sang = sin(ang)
        1447           0 :   cang = cos(ang)
        1448           0 :   cte = ctg*std + stg*ctd*cb
        1449           0 :   ste = sqrt(1._r8 - cte*cte)
        1450           0 :   sa = sb*ctd/ste
        1451           0 :   ca = (std*stg - ctd*ctg*cb)/ste
        1452           0 :   capang = ca*cang - sa*sang
        1453           0 :   sapang = ca*sang + sa*cang
        1454           0 :   stfcpa = ste*ctp*capang - cte*stp
        1455           0 :   stfspa = sapang*ste
        1456           0 :   alon = atan2(stfspa,stfcpa)*rtd
        1457           0 :   r = alt + re
        1458           0 :   ha = alt + r*cottd*cottd
        1459           0 :   a = 1._r8 + ha/req
        1460             :
        1461             : end subroutine dipapx
        1462             : !-----------------------------------------------------------------------
        1463           0 : subroutine itrace(iapx)
        1464             :   save
        1465             : !
        1466             : ! Uses 4-point ADAMS formula after initialization.
        1467             : ! First 7 iterations advance point by 3 steps.
        1468             : !
        1469             : ! y(3), yp(3), yapx(3,3), sgn and nstp are in module data
        1470             : ! yploc(3,4) is local
        1471             : !
        1472             : ! Arg:
        1473             :   integer,intent(out) :: iapx
        1474             : !
        1475             : ! Local:
        1476             :   integer :: i,j
        1477             :   real(r8) :: yploc(3,4) ! local yp (i.e., not module data yp)
        1478             :   real(r8) :: term,d2,d6,d12,d24,rc,rp
        1479             :
        1480           0 :   iapx = 1
        1481             : !
        1482             : ! Field line is defined by the following differential equations
        1483             : ! in cartesian coordinates.
        1484             : ! (yapx,yp,y are module data)
        1485             : !
        1486           0 :   yploc(1,4) = sgn*bx/bb
        1487           0 :   yploc(2,4) = sgn*by/bb
        1488           0 :   yploc(3,4) = sgn*bz/bb
        1489             :
        1490           0 :   if (nstp > 7) then
        1491           0 :     do i=1,3
        1492           0 :       yapx(i,1) = yapx(i,2)
        1493           0 :       yapx(i,2) = y(i)
        1494           0 :       yp(i) = y(i)
        1495           0 :       term = 55._r8*yploc(i,4)-59._r8*yploc(i,3)+37._r8*yploc(i,2)-9._r8*yploc(i,1)
        1496           0 :       y(i) = yp(i) + d24*term
        1497           0 :       yapx(i,3) = y(i)
        1498           0 :       do j=1,3
        1499           0 :         yploc(i,j) = yploc(i,j+1)
        1500             :       enddo
        1501             :     enddo
        1502           0 :     rc = rdus ( y(1),  y(2),  y(3))
        1503           0 :     rp = rdus (yp(1), yp(2), yp(3))
        1504           0 :     if (rc < rp) iapx=2
        1505           0 :     return
        1506             :   endif
        1507             :
        1508           0 :   do i=1,3
        1509           0 :     select case (nstp)
        1510             :       case (1)
        1511           0 :         d2  = ds/2._r8
        1512           0 :         d6  = ds/6._r8
        1513           0 :         d12 = ds/12._r8
        1514           0 :         d24 = ds/24._r8
        1515           0 :         yploc(i,1)= yploc(i,4)
        1516           0 :         yp(i)     = y(i)
        1517           0 :         yapx(i,1) = y(i)
        1518           0 :         y(i) = yp(i) + ds*yploc(i,1)
        1519             :       case (2)
        1520           0 :         yploc(i,2) = yploc(i,4)
        1521           0 :         y(i) = yp(i) + d2*(yploc(i,2)+yploc(i,1))
        1522             :       case (3)
        1523           0 :         y(i) = yp(i) + d6*(2._r8*yploc(i,4)+yploc(i,2)+3._r8*yploc(i,1))
        1524             :       case (4)
        1525           0 :         yploc(i,2) = yploc(i,4)
        1526           0 :         yapx(i,2)  = y(i)
        1527           0 :         yp(i)      = y(i)
        1528           0 :         y(i)       = yp(i) + d2*(3._r8*yploc(i,2)-yploc(i,1))
        1529             :       case (5)
        1530           0 :         y(i) = yp(i) + d12*(5._r8*yploc(i,4)+8._r8*yploc(i,2)-yploc(i,1))
        1531             :       case (6)
        1532           0 :         yploc(i,3) = yploc(i,4)
        1533           0 :         yp(i)      = y(i)
        1534           0 :         yapx(i,3)  = y(i)
        1535           0 :         y(i)       = yp(i) + d12*(23._r8*yploc(i,3)-16._r8*yploc(i,2)+5._r8*yploc(i,1))
        1536             :       case (7)
        1537           0 :         yapx(i,1) = yapx(i,2)
        1538           0 :         yapx(i,2) = yapx(i,3)
        1539           0 :         y(i) = yp(i) + d24*(9._r8*yploc(i,4)+19._r8*yploc(i,3)-5._r8*yploc(i,2)+yploc(i,1))
        1540           0 :         yapx(i,3) = y(i)
        1541             :       case default
        1542           0 :         write(iulog,"('>>> itrace: unresolved case nstp=',i4)") nstp
        1543           0 :         call endrun( 'itrace' )
        1544             :     end select
        1545             :   enddo
        1546             : !
        1547             : ! Signal if apex passed:
        1548             : !
        1549           0 :   if (nstp == 6 .or. nstp == 7) then
        1550           0 :     rc = rdus( yapx(1,3), yapx(2,3), yapx(3,3))
        1551           0 :     rp = rdus( yapx(1,2), yapx(2,2), yapx(3,2))
        1552           0 :     if (rc < rp) iapx=2
        1553             :   endif
        1554             :
        1555             : end subroutine itrace
        1556             : !-----------------------------------------------------------------------
        1557           0 : real(r8) function rdus(d,e,f)
        1558             :   real(r8),intent(in) :: d,e,f
        1559           0 :   rdus = sqrt(d**2 + e**2 + f**2)
        1560           0 : end function rdus
        1561             : !-----------------------------------------------------------------------
        1562           0 : subroutine fndapx(alt,zmag,a,alat,alon)
        1563             : !
        1564             : ! Find apex coords once tracing has signalled that the apex has been passed.
        1565             : !
        1566             : ! Args:
        1567             :   real(r8),intent(in) :: alt,zmag
        1568             :   real(r8),intent(out) :: a,alat,alon
        1569             : !
        1570             : ! Local:
        1571             :   integer :: i,iflag_convrt, iflag_feldg
        1572             :   real(r8) :: z(3),ht(3),yloc(3),gdlt,gdln,x,ydum,f,rho,xinter,rasq,xlon,ang,&
        1573             :     cang,sang,r,cte,ste,stfcpa,stfspa
        1574             : !
        1575             : ! Get geodetic field components.
        1576             : !
        1577           0 :   iflag_feldg = 1
        1578           0 :   iflag_convrt = 3
        1579           0 :   do i=1,3
        1580           0 :     rho = sqrt(yapx(1,i)**2+yapx(2,i)**2)
        1581           0 :     call convrt(iflag_convrt,gdlt,ht(i),rho,yapx(3,i))
        1582           0 :     gdln = rtd*atan2(yapx(2,i),yapx(1,i))
        1583           0 :     call feldg(iflag_feldg,gdlt,gdln,ht(i),x,ydum,z(i),f)
        1584             :   enddo
        1585             : !
        1586             : ! Find cartesian coordinates at dip equator by interpolation
        1587             : !
        1588           0 :   do i=1,3
        1589           0 :     call fint(z(1),z(2),z(3),yapx(i,1),yapx(i,2),yapx(i,3),0._r8,yloc(i))
        1590             :   enddo
        1591             : !
        1592             : ! Find apex height by interpolation
        1593             : !
        1594           0 :   call fint(z(1),z(2),z(3),ht(1),ht(2),ht(3),0._r8,xinter)
        1595             : !
        1596             : ! Ensure that XINTER is not less than original starting altitude:
        1597           0 :   xinter = max(alt,xinter)
        1598           0 :   a = (req+xinter)/req
        1599             : !
        1600             : ! Find apex coordinates , giving alat sign of dip at starting point.
        1601             : ! Alon is the value of the geomagnetic longitude at the apex.
        1602             : !
        1603           0 :   if (a < 1._r8) then
        1604           0 :     write(iulog,"('>>> fndapx: a=',e12.4,' < 1.')") a
        1605           0 :     call endrun( 'fndapx' )
        1606             :   endif
        1607             :
        1608           0 :   rasq = rtd*acos(sqrt(1._r8/a))
        1609           0 :   alat = sign(rasq,zmag)
        1610             : !
        1611             : ! Algorithm for ALON:
        1612             : !   Use spherical coordinates.
        1613             : !   Let GP be geographic pole.
        1614             : !   Let GM be geomagnetic pole (colatitude COLAT, east longitude ELON).
        1615             : !   Let XLON be longitude of apex.
        1616             : !   Let TE be colatitude of apex.
        1617             : !   Let ANG be longitude angle from GM to apex.
        1618             : !   Let TP be colatitude of GM.
        1619             : !   Let TF be arc length between GM and apex.
        1620             : !   Let PA = ALON be geomagnetic longitude, i.e., Pi minus angle measured
        1621             : !     counterclockwise from arc GM-apex to arc GM-GP.
        1622             : !   Then, using notation C=cos, S=sin, spherical-trigonometry formulas
        1623             : !     for the functions of the angles are as shown below.  Note: STFCPA,
        1624             : !     STFSPA are sin(TF) times cos(PA), sin(PA), respectively.
        1625             : !
        1626           0 :   xlon = atan2(yloc(2),yloc(1))
        1627           0 :   ang  = xlon-elon*dtr
        1628           0 :   cang = cos(ang)
        1629           0 :   sang = sin(ang)
        1630           0 :   r    = sqrt(yloc(1)**2+yloc(2)**2+yloc(3)**2)
        1631           0 :   cte  = yloc(3)/r
        1632           0 :   ste  = sqrt(1._r8-cte*cte)
        1633           0 :   stfcpa = ste*ctp*cang - cte*stp
        1634           0 :   stfspa = sang*ste
        1635           0 :   alon = atan2(stfspa,stfcpa)*rtd
        1636             :
        1637           0 : end subroutine fndapx
        1638             : !-----------------------------------------------------------------------
        1639           0 : subroutine fint(a1,a2,a3,a4,a5,a6,a7,result)
        1640             : !
        1641             : ! Second degree interpolation
        1642             : !
        1643             : ! Args:
        1644             :   real(r8),intent(in) :: a1,a2,a3,a4,a5,a6,a7
        1645             :   real(r8),intent(out) :: result
        1646             :
        1647             :   result = ((a2-a3)*(a7-a2)*(a7-a3)*a4-(a1-a3)*(a7-a1)*(a7-a3)*a5+ &
        1648           0 :     (a1-a2)*(a7-a1)*(a7-a2)*a6)/((a1-a2)*(a1-a3)*(a2-a3))
        1649             :
        1650           0 : end subroutine fint
        1651             : !-----------------------------------------------------------------------
        1652           0 : subroutine gm2gc(gmlat,gmlon,gclat,gclon)
        1653             : !
        1654             : ! Args:
        1655             :   real(r8),intent(in)  :: gmlat,gmlon
        1656             :   real(r8),intent(out) :: gclat,gclon
        1657             : !
        1658             : ! Local:
        1659             :   real(r8) :: stm,ctm,ctc
        1660             :
        1661           0 :   stm = cos(gmlat*dtr)
        1662           0 :   ctm = sin(gmlat*dtr)
        1663           0 :   ctc = ctp*ctm - stp*stm*cos(gmlon*dtr) ! ctp,stp are module data
        1664           0 :   ctc = min(ctc,1._r8)
        1665           0 :   ctc = max(ctc,-1._r8)
        1666           0 :   gclat = asin(ctc)*rtd
        1667           0 :   gclon = atan2(stp*stm*sin(gmlon*dtr),ctm-ctp*ctc)
        1668             : !
        1669             : ! elon is in module data, and was set by dypol (called from apex_mka)
        1670             : !
        1671           0 :   gclon = gclon*rtd + elon
        1672           0 :   if (gclon < -180._r8) gclon = gclon + 360._r8
        1673             :
        1674           0 : end subroutine gm2gc
        1675             : !-----------------------------------------------------------------------
        1676           0 : subroutine intrp(glat,glon,alt, gplat,gplon,gpalt, nlat,nlon,nalt, &
        1677             :                  fx,fy,fz,fv,                                      &
        1678             :                  dfxdth,dfydth,dfzdth,dfvdth,                      &
        1679             :                  dfxdln,dfydln,dfzdln,dfvdln,                      &
        1680             :                  dfxdh ,dfydh ,dfzdh ,dfvdh, ier)
        1681             : !
        1682             : ! Args:
        1683             : !
        1684             :   real(r8),intent(in)    :: glat,glon,alt
        1685             :   integer,intent(in) :: nlat,nlon,nalt
        1686             :   real(r8),intent(in)    :: gplat(nlat),gplon(nlon),gpalt(nalt)
        1687             :   real(r8),intent(out)   ::          &
        1688             :     fx,fy,fz,fv,                 &
        1689             :     dfxdth,dfydth,dfzdth,dfvdth, &
        1690             :     dfxdln,dfydln,dfzdln,dfvdln, &
        1691             :     dfxdh ,dfydh ,dfzdh ,dfvdh
        1692             :   integer,intent(out) :: ier
        1693             : !
        1694             : ! Local:
        1695             : !
        1696             :   integer :: i,j,k,i0,j0,k0
        1697             :   real(r8) :: glonloc,xi,dlon,yj,dlat,hti,diht,zk,fac,omfac
        1698             :   real(r8) :: dfxdn,dfxde,dfxdd, &
        1699             :           dfydn,dfyde,dfydd, &
        1700             :           dfzdn,dfzde,dfzdd, &
        1701             :           dfvdn,dfvde,dfvdd, &
        1702             :           dmf,dmdfdn,dmdfde,dmdfdd
        1703             :
        1704           0 :   ier = 0
        1705           0 :   glonloc = glon
        1706           0 :   if (glonloc < gplon(1))    glonloc = glonloc + 360._r8
        1707           0 :   if (glonloc > gplon(nlon)) glonloc = glonloc - 360._r8
        1708             : !
        1709           0 :   i0 = 0
        1710           0 :   do i=1,nlat-1
        1711           0 :     if (glat >= gplat(i).and.glat <= gplat(i+1)) then
        1712           0 :       i0 = i
        1713           0 :       dlat = gplat(i+1)-gplat(i)
        1714           0 :       xi = (glat - gplat(i)) / dlat
        1715           0 :       exit
        1716             :     endif
        1717             :   enddo
        1718           0 :   if (i0==0) then
        1719             :     write(iulog,"('>>> intrp: could not bracket glat=',f9.3,' in gplat=',/,(6f9.2))") &
        1720           0 :       glat,gplat
        1721           0 :     ier = 1
        1722           0 :     return
        1723             :   endif
        1724             :
        1725           0 :   j0 = 0
        1726           0 :   do j=1,nlon-1
        1727           0 :     if (glon >= gplon(j).and.glon <= gplon(j+1)) then
        1728           0 :       j0 = j
        1729           0 :       dlon = gplon(j+1)-gplon(j)
        1730           0 :       yj = (glon - gplon(j)) / dlon
        1731           0 :       exit
        1732             :     endif
        1733             :   enddo
        1734           0 :   if (j0==0) then
        1735             :     write(iulog,"('>>> intrp: could not bracket glon=',f9.3,' in gplon=',/,(6f9.2))") &
        1736           0 :       glon,gplon
        1737           0 :     ier = 1
        1738           0 :     return
        1739             :   endif
        1740             :
        1741           0 :   k0 = 0
        1742           0 :   do k=1,nalt-1
        1743           0 :     if (alt >= gpalt(k).and.alt <= gpalt(k+1)) then
        1744           0 :       k0 = k
        1745           0 :       hti = re/(re+alt)
        1746           0 :       diht = re/(re+gpalt(k+1)) - re/(re+gpalt(k))
        1747           0 :       zk = (hti - re/(re+gpalt(k))) / diht
        1748           0 :       exit
        1749             :     endif
        1750             :   enddo
        1751           0 :   if (k0==0) then
        1752             :     write(iulog,"('>>> intrp: could not bracket alt=',f12.3,' in gpalt=',/,(6f12.2))") &
        1753           0 :       alt,gpalt
        1754           0 :     ier = 1
        1755           0 :     return
        1756             :   endif
        1757             :
        1758           0 :   call trilin(xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fx,dfxdn,dfxde,dfxdd)
        1759           0 :   dfxdth = -dfxdn*rtd/dlat
        1760           0 :   dfxdln =  dfxde*rtd/dlon
        1761           0 :   dfxdh  = -hti*hti*dfxdd/(re*diht)
        1762             :
        1763           0 :   call trilin(yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fy,dfydn,dfyde,dfydd)
        1764           0 :   dfydth = -dfydn*rtd/dlat
        1765           0 :   dfydln =  dfyde*rtd/dlon
        1766           0 :   dfydh  = -hti*hti*dfydd/(re*diht)
        1767             :
        1768           0 :   call trilin(zarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fz,dfzdn,dfzde,dfzdd)
        1769           0 :   dfzdth = -dfzdn*rtd/dlat
        1770           0 :   dfzdln =  dfzde*rtd/dlon
        1771           0 :   dfzdh  = -hti*hti*dfzdd/(re*diht)
        1772             :
        1773           0 :   call trilin(varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fv,dfvdn,dfvde,dfvdd)
        1774           0 :   dfvdth = -dfvdn*rtd/dlat
        1775           0 :   dfvdln =  dfvde*rtd/dlon
        1776           0 :   dfvdh  = -hti*hti*dfvdd/(re*diht)
        1777             :
        1778           0 :   if (nlat < 3) return
        1779             : !
        1780             : ! Improve calculation of longitudinal derivatives near poles
        1781             : !
        1782           0 :   if (glat < dlat-90._r8) then
        1783           0 :     fac = .5_r8*xi
        1784           0 :     omfac = 1._r8 - fac
        1785           0 :     xi = xi - 1._r8
        1786           0 :     i0 = i0 + 1
        1787           0 :     call trilin (xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1788           0 :     dfxdln = dfxdln*omfac + fac*dmdfde*rtd/dlon
        1789           0 :     call trilin (yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1790           0 :     dfydln = dfydln*omfac + fac*dmdfde*rtd/dlon
        1791           0 :     call trilin (varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1792           0 :     dfvdln = dfvdln*omfac + fac*dmdfde*rtd/dlon
        1793             :   endif
        1794             :
        1795           0 :   if (glat > 90._r8-dlat) then
        1796           0 :     fac = .5_r8*(1._r8-xi)
        1797           0 :     omfac = 1._r8 - fac
        1798           0 :     xi = xi + 1._r8
        1799           0 :     i0 = i0 - 1
        1800           0 :     call trilin (xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1801           0 :     dfxdln = dfxdln*omfac + fac*dmdfde*rtd/dlon
        1802           0 :     call trilin (yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1803           0 :     dfydln = dfydln*omfac + fac*dmdfde*rtd/dlon
        1804           0 :     call trilin (varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
        1805           0 :     dfvdln = dfvdln*omfac + fac*dmdfde*rtd/dlon
        1806             :   endif
        1807             :
        1808             : end subroutine intrp
        1809             : !-----------------------------------------------------------------------
        1810           0 : subroutine trilin(u,xi,yj,zk,fu,dfudx,dfudy,dfudz)
        1811             : !
        1812             : ! Args:
        1813             :   real(r8),intent(in) :: &
        1814             :     u(1:2,1:2,1:2),      & ! u(1,1,1) is address of lower corner of interpolation box
        1815             :     xi,  & ! fractional distance across box in x direction
        1816             :     yj,  & ! fractional distance across box in y direction
        1817             :     zk     ! fractional distance across box in z direction
        1818             :   real(r8),intent(out)   :: &
        1819             :     fu,                 & ! interpolated value of u
        1820             :     dfudx,              & ! interpolated derivative of u with respect to i (x direction)
        1821             :     dfudy,              & ! interpolated derivative of u with respect to j (y direction)
        1822             :     dfudz                 ! interpolated derivative of u with respect to k (z direction)
        1823             : !
        1824             : ! Local:
        1825             :   real(r8) :: omxi,omyj,omzk
        1826             :
        1827             : ! write(iulog,"('Enter trilin: xi,yj,zk=',3e12.4)") xi,yj,zk
        1828             : ! write(iulog,"('Enter trilin: u(1,1,1),u(1,2,1),u(1,1,2),u(1,2,2)=',4e12.4)") &
        1829             : !                          u(1,1,1),u(1,2,1),u(1,1,2),u(1,2,2)
        1830             : ! write(iulog,"('Enter trilin: u(2,1,1),u(2,2,1),u(2,1,2),u(2,2,2)=',4e12.4)") &
        1831             : !                          u(2,1,1),u(2,2,1),u(2,1,2),u(2,2,2)
        1832             :
        1833           0 :   omxi = 1._r8 - xi
        1834           0 :   omyj = 1._r8 - yj
        1835           0 :   omzk = 1._r8 - zk
        1836             :
        1837             :   fu = u(1,1,1)*omxi*omyj*omzk &
        1838             :      + u(2,1,1)*xi*omyj*omzk   &
        1839             :      + u(1,2,1)*omxi*yj*omzk   &
        1840             :      + u(1,1,2)*omxi*omyj*zk   &
        1841             :      + u(2,2,1)*xi*yj*omzk     &
        1842             :      + u(2,1,2)*xi*omyj*zk     &
        1843             :      + u(1,2,2)*omxi*yj*zk     &
        1844           0 :      + u(2,2,2)*xi*yj*zk
        1845             :
        1846             :   dfudx = (u(2,1,1)-u(1,1,1))*omyj*omzk &
        1847             :         + (u(2,2,1)-u(1,2,1))*yj*omzk   &
        1848             :         + (u(2,1,2)-u(1,1,2))*omyj*zk   &
        1849           0 :         + (u(2,2,2)-u(1,2,2))*yj*zk
        1850             :   dfudy = (u(1,2,1)-u(1,1,1))*omxi*omzk &
        1851             :         + (u(2,2,1)-u(2,1,1))*xi*omzk   &
        1852             :         + (u(1,2,2)-u(1,1,2))*omxi*zk   &
        1853           0 :         + (u(2,2,2)-u(2,1,2))*xi*zk
        1854             :   dfudz = (u(1,1,2)-u(1,1,1))*omxi*omyj &
        1855             :         + (u(2,1,2)-u(2,1,1))*xi*omyj   &
        1856             :         + (u(1,2,2)-u(1,2,1))*omxi*yj   &
        1857           0 :         + (u(2,2,2)-u(2,2,1))*xi*yj
        1858             :
        1859           0 : end subroutine trilin
        1860             : !-----------------------------------------------------------------------
        1861           0 : subroutine adpl(glat,glon,cth,sth,fx,fy,fz,fv, &
        1862             :                 dfxdth,dfydth,dfzdth,dfvdth,dfxdln,dfydln,dfzdln,dfvdln)
        1863             : !
        1864             : !  Add-back of pseudodipole component to x,y,z,v and their derivatives.
        1865             : !
        1866             : ! Args:
        1867             :   real(r8),intent(in)      :: glat,glon
        1868             :   real(r8),intent(out)     :: cth,sth
        1869             :   real(r8),intent(inout)   ::    &
        1870             :      fx,fy,fz,fv,                &
        1871             :     dfxdth,dfydth,dfzdth,dfvdth, &
        1872             :     dfxdln,dfydln,dfzdln,dfvdln
        1873             : !
        1874             : ! Local:
        1875             :   real(r8) :: cph,sph,ctm
        1876             :
        1877           0 :   cph = cos((glon-elon)*dtr)
        1878           0 :   sph = sin((glon-elon)*dtr)
        1879           0 :   cth = sin(glat*dtr)
        1880           0 :   sth = cos(glat*dtr)
        1881           0 :   ctm = ctp*cth + stp*sth*cph
        1882           0 :   fx = fx + sth*ctp*cph - cth*stp
        1883           0 :   fy = fy + sth*sph
        1884           0 :   fz = fz + ctm
        1885           0 :   fv = fv - ctm
        1886             :
        1887           0 :   dfxdth = dfxdth + ctp*cth*cph + stp*sth
        1888           0 :   dfydth = dfydth + cth*sph
        1889           0 :   dfzdth = dfzdth - ctp*sth + stp*cth*cph
        1890           0 :   dfvdth = dfvdth + ctp*sth - stp*cth*cph
        1891             :
        1892           0 :   dfxdln = dfxdln - ctp*sth*sph
        1893           0 :   dfydln = dfydln + sth*cph
        1894           0 :   dfzdln = dfzdln - stp*sth*sph
        1895           0 :   dfvdln = dfvdln + stp*sth*sph
        1896             :
        1897           0 : end subroutine adpl
        1898             : !-----------------------------------------------------------------------
        1899           0 : subroutine setmiss(xmiss,xlatm,alon,vmp,b,bmag,be3,sim,si,f,d,w, &
        1900             :   bhat,d1,d2,d3,e1,e2,e3,f1,f2)
        1901             : !
        1902             : ! Args:
        1903             :   real(r8),intent(in)  :: xmiss
        1904             :   real(r8),intent(out) :: xlatm,alon,vmp,bmag,be3,sim,si,f,d,w
        1905             :   real(r8),dimension(3),intent(out) :: bhat,d1,d2,d3,e1,e2,e3,b
        1906             :   real(r8),dimension(2),intent(out) :: f1,f2
        1907             :
        1908           0 :   xlatm = xmiss
        1909           0 :   alon  = xmiss
        1910           0 :   vmp   = xmiss
        1911           0 :   bmag  = xmiss
        1912           0 :   be3   = xmiss
        1913           0 :   sim   = xmiss
        1914           0 :   si    = xmiss
        1915           0 :   f     = xmiss
        1916           0 :   d     = xmiss
        1917           0 :   w     = xmiss
        1918           0 :   bhat  = xmiss
        1919           0 :   d1    = xmiss
        1920           0 :   d2    = xmiss
        1921           0 :   d3    = xmiss
        1922           0 :   e1    = xmiss
        1923           0 :   e2    = xmiss
        1924           0 :   e3    = xmiss
        1925           0 :   b     = xmiss
        1926           0 :   f1    = xmiss
        1927           0 :   f2    = xmiss
        1928             :
        1929           0 : end subroutine setmiss
        1930             : !-----------------------------------------------------------------------
        1931           0 : subroutine cofrm(date)
        1932             :   implicit none
        1933             : !
        1934             : ! Input arg:
        1935             :   real(r8),intent(in) :: date
        1936             : !
        1937             : ! Local:
        1938             :   integer :: m,n,i,l,ll,lm,nmx,nc,kmx,k,nn
        1939             :   real(r8) :: t,one,tc,f,f0
        1940             :
        1941             :   integer :: ngh !  = n1*ncn1 + n2*ncn2 + 1
        1942           0 :   real(r8) :: gh(n1*ncn1 + n2*ncn2 + 1)
        1943             :
        1944             :   real(r8),parameter :: alt = 0._r8
        1945             :   integer, parameter :: isv=0
        1946             :
        1947           0 :   ngh = n1*ncn1 + n2*ncn2 + 1 ! not sure why the extra +1
        1948             :
        1949           0 :   if (date < apex_beg_yr .or. date > apex_end_yr) then
        1950           0 :     write(iulog,"('>>> cofrm: date=',f8.2,' Date must be >= ',I4,' and <= ',I4)") date,apex_beg_yr,apex_end_yr+5
        1951           0 :     call endrun( 'cofrm' )
        1952             :   endif
        1953           0 :   if (date > apex_end_yr-5) then
        1954           0 :    if (masterproc .and. .not. first_warning) then
        1955           0 :     write(iulog,"('>>> WARNING cofrm:')")
        1956           0 :     write(iulog,"(/,'   This version of IGRF is intended for use up to ')")
        1957           0 :     write(iulog,"('     2020. Values for ',f9.3,' will be computed but')") date
        1958           0 :     write(iulog,"('     may be of reduced accuracy.',/)")
        1959           0 :     first_warning=.true.
        1960             :    endif
        1961             :   endif
        1962             : !
        1963             : ! Set gh from g1,g2:
        1964             : !
        1965           0 :   do n=1,ncn1
        1966           0 :     i = (n-1)*n1
        1967           0 :     gh(i+1:i+n1) = g1(:,n)
        1968             : !   write(iulog,"('cofrm: n=',i3,' i+1:i+n1=',i4,':',i4)") n,i+1,i+n1
        1969             :   enddo
        1970           0 :   do n=1,ncn2
        1971           0 :     i = n1*ncn1 + (n-1)*n2
        1972           0 :     gh(i+1:i+n2) = g2(:,n)
        1973             : !   write(iulog,"('cofrm: n=',i3,' i+1:i+n2=',i4,':',i4)") n,i+1,i+n2
        1974             :   enddo
        1975           0 :   gh(ngh) = 0._r8 ! not sure why gh is dimensioned with the extra element, so set it to 0.
        1976             :
        1977           0 :   if (date < apex_end_yr-10) then
        1978           0 :     t   = 0.2_r8*(date - year1)
        1979           0 :     ll  = t
        1980           0 :     one = ll
        1981           0 :     t   = t - one
        1982           0 :     if (date < year2-5) then
        1983           0 :       nmx   = 10
        1984           0 :       nc    = nmx*(nmx+2)
        1985           0 :       ll    = nc*ll
        1986           0 :       kmx   = (nmx+1)*(nmx+2)/2
        1987             :     else
        1988           0 :       nmx   = 13
        1989           0 :       nc    = nmx*(nmx+2)
        1990           0 :       ll    = 0.2_r8*(date - (year2-5))
        1991           0 :       ll    = 120*19 + nc*ll
        1992           0 :       kmx   = (nmx+1)*(nmx+2)/2
        1993             :     endif
        1994           0 :     tc    = 1.0_r8 - t
        1995             :     if (isv.eq.1) then
        1996             :       tc = -0.2_r8
        1997             :       t = 0.2_r8
        1998             :     endif
        1999             :   else ! date >= apex_end_yr-10
        2000           0 :     t     = date - (apex_end_yr-10)
        2001           0 :     tc    = 1.0_r8
        2002             :     if (isv.eq.1) then
        2003             :       t = 1.0_r8
        2004             :       tc = 0.0_r8
        2005             :     end if
        2006           0 :     ll    = n1*ncn1 + n2*(ncn2-2) ! corresponds to apex_end_yr-10
        2007           0 :     nmx   = 13
        2008           0 :     nc    = nmx*(nmx+2)
        2009           0 :     kmx   = (nmx+1)*(nmx+2)/2
        2010             :   endif ! date < apex_end_yr-10
        2011           0 :   l = 1
        2012           0 :   m = 1
        2013           0 :   n = 0
        2014             : !
        2015             : ! Set outputs gb(ncoef) and gv(ncoef)
        2016             : ! These are module data above.
        2017             : !
        2018           0 :   gb(:) = 0._r8
        2019           0 :   gv(:) = 0._r8
        2020           0 :   f0 = -1.e-5_r8
        2021           0 :   do k=2,kmx
        2022           0 :     if (n < m) then
        2023           0 :       m = 0
        2024           0 :       n = n+1
        2025             :     endif ! n < m
        2026           0 :     lm = ll + l
        2027           0 :     if (m == 0) f0 = f0 * dble(n)/2._r8
        2028           0 :     if (m == 0) f  = f0 / sqrt(2.0_r8)
        2029           0 :     nn = n+1
        2030             :
        2031           0 :     if (m /= 0) then
        2032           0 :       f = f / sqrt(dble(n-m+1) / dble(n+m) )
        2033           0 :       gb(l+1)  = (tc*gh(lm) + t*gh(lm+nc))* f
        2034             :     else
        2035           0 :       gb(l+1)  = (tc*gh(lm) + t*gh(lm+nc))* f0
        2036             :     endif
        2037           0 :     gv(l+1) = gb(l+1)/dble(nn)
        2038           0 :     if (m /= 0) then
        2039           0 :       gb(l+2)  = (tc*gh(lm+1) + t*gh(lm+nc+1))*f
        2040           0 :       gv(l+2) = gb(l+2)/dble(nn)
        2041           0 :       l = l+2
        2042             :     else
        2043             :       l = l+1
        2044             :     endif
        2045           0 :     m = m+1
        2046             :   enddo
        2047             :
        2048             : ! write(iulog,"('cofrm: ncoef=',i4,' gb=',/,(6f12.3))") ncoef,gb
        2049             : ! write(iulog,"('cofrm: ncoef=',i4,' gv=',/,(6f12.3))") ncoef,gv
        2050             :
        2051           0 : end subroutine cofrm
        2052             : !-----------------------------------------------------------------------
        2053           0 : subroutine apex_subsol(iyr,iday,ihr,imn,sec,sbsllat,sbsllon)
        2054             : !
        2055             : ! Find subsolar geographic latitude and longitude given the
        2056             : ! date and time (Universal Time).
        2057             : !
        2058             : ! This is based on formulas in Astronomical Almanac for the
        2059             : ! year 1996, p.  C24. (U.S.  Government Printing Office,
        2060             : ! 1994).  According to the Almanac, results are good to at
        2061             : ! least 0.01 degree latitude and 0.025 degree longitude
        2062             : ! between years 1950 and 2050.  Accuracy for other years has
        2063             : ! not been tested although the algorithm has been designed to
        2064             : ! accept input dates from 1601 to 2100.  Every day is assumed
        2065             : ! to have exactly 86400 seconds; thus leap seconds that
        2066             : ! sometimes occur on June 30 and December 31 are ignored:
        2067             : ! their effect is below the accuracy threshold of the algorithm.
        2068             : !
        2069             : ! 961026 A. D. Richmond, NCAR
        2070             : !
        2071             : ! Input Args:
        2072             :   integer,intent(in) :: &
        2073             :     iyr,   & ! Year (e.g., 1994). IYR must be in the range: 1601 to 2100.
        2074             :     iday,  & ! Day number of year (e.g., IDAY = 32 for Feb 1)
        2075             :     ihr,   & ! Hour of day    (e.g., 13 for 13:49)
        2076             :     imn      ! Minute of hour (e.g., 49 for 13:49)
        2077             :   real(r8),intent(in) :: sec ! Second and fraction after the hour/minute.
        2078             : !
        2079             : ! Output Args:
        2080             :   real(r8),intent(out) :: &
        2081             :     sbsllat, & ! geographic latitude of subsolar point (degrees)
        2082             :     sbsllon    ! geographic longitude of subsolar point (-180 to +180)
        2083             : !
        2084             : ! Local:
        2085             :   integer,parameter :: minyear=1601, maxyear = 2100
        2086             :   real(r8),parameter :: & ! Use local params for compatability w/ legacy code,
        2087             :                           ! but probably would be ok to use module data dtr,rtd
        2088             :     d2r=0.0174532925199432957692369076847_r8, &
        2089             :     r2d=57.2957795130823208767981548147_r8
        2090             :   real(r8) :: yr,l0,g0,ut,df,lf,gf,l,g,grad,n,epsilon,epsrad,alpha,delta,&
        2091             :     etdeg,aptime,lambda,lamrad,sinlam
        2092             :   integer :: nleap,ncent,nrot
        2093             :
        2094           0 :   sbsllat=0._r8 ; sbsllon=0._r8
        2095             :
        2096           0 :   yr = iyr-2000
        2097             : !
        2098             : ! nleap (final) = number of leap days from (2000 January 1) to (IYR January 1)
        2099             : !                 (negative if iyr is before 1997)
        2100           0 :   nleap = (iyr-1601)/4
        2101           0 :   nleap = nleap - 99
        2102           0 :   if (iyr <= year1) then
        2103           0 :     if (iyr < minyear) then
        2104           0 :       write(iulog,*) 'subsolr invalid before ',minyear,': input year = ',iyr
        2105           0 :       call endrun( 'subsolr' )
        2106             :     endif
        2107           0 :     ncent = (iyr-minyear)/100
        2108           0 :     ncent = 3 - ncent
        2109           0 :     nleap = nleap + ncent
        2110             :   endif
        2111           0 :   if (iyr > maxyear) then
        2112           0 :     write(iulog,*) 'subsolr invalid after ',maxyear,':  input year = ',iyr
        2113           0 :     call endrun( 'subsolr' )
        2114             :   endif
        2115             : !
        2116             : ! L0 = Mean longitude of Sun at 12 UT on January 1 of IYR:
        2117             : !     L0 = 280.461 + .9856474*(365*(YR-NLEAP) + 366*NLEAP)
        2118             : !          - (ARBITRARY INTEGER)*360.
        2119             : !        = 280.461 + .9856474*(365*(YR-4*NLEAP) + (366+365*3)*NLEAP)
        2120             : !          - (ARBITRARY INTEGER)*360.
        2121             : !        = (280.461 - 360.) + (.9856474*365 - 360.)*(YR-4*NLEAP)
        2122             : !          + (.9856474*(366+365*3) - 4*360.)*NLEAP,
        2123             : !  where ARBITRARY INTEGER = YR+1.  This gives:
        2124             : !
        2125           0 :       l0 = -79.549_r8 + (-.238699_r8*(yr-4*nleap) + 3.08514e-2_r8*nleap)
        2126             : !
        2127             : ! G0 = Mean anomaly at 12 UT on January 1 of IYR:
        2128             : !     G0 = 357.528 + .9856003*(365*(YR-NLEAP) + 366*NLEAP)
        2129             : !          - (ARBITRARY INTEGER)*360.
        2130             : !        = 357.528 + .9856003*(365*(YR-4*NLEAP) + (366+365*3)*NLEAP)
        2131             : !          - (ARBITRARY INTEGER)*360.
        2132             : !        = (357.528 - 360.) + (.9856003*365 - 360.)*(YR-4*NLEAP)
        2133             : !          + (.9856003*(366+365*3) - 4*360.)*NLEAP,
        2134             : !  where ARBITRARY INTEGER = YR+1.  This gives:
        2135             : !
        2136           0 :       g0 = -2.472_r8 + (-.2558905_r8*(yr-4*nleap) - 3.79617e-2_r8*nleap)
        2137             : !
        2138             : ! Universal time in seconds:
        2139           0 :       ut = dble(ihr*3600 + imn*60) + sec
        2140             : !
        2141             : ! Days (including fraction) since 12 UT on January 1 of IYR:
        2142           0 :       df = (ut/86400._r8 - 1.5_r8) + iday
        2143             : !
        2144             : ! Addition to Mean longitude of Sun since January 1 of IYR:
        2145           0 :       lf = .9856474_r8*df
        2146             : !
        2147             : ! Addition to Mean anomaly since January 1 of IYR:
        2148           0 :       gf = .9856003_r8*df
        2149             : !
        2150             : ! Mean longitude of Sun:
        2151           0 :       l = l0 + lf
        2152             : !
        2153             : ! Mean anomaly:
        2154           0 :       g = g0 + gf
        2155           0 :       grad = g*d2r
        2156             : !
        2157             : ! Ecliptic longitude:
        2158           0 :       lambda = l + 1.915_r8*sin(grad) + .020_r8*sin(2._r8*grad)
        2159           0 :       lamrad = lambda*d2r
        2160           0 :       sinlam = sin(lamrad)
        2161             : !
        2162             : ! Days (including fraction) since 12 UT on January 1 of 2000:
        2163           0 :       n = df + 365._r8*yr + dble(nleap)
        2164             : !
        2165             : ! Obliquity of ecliptic:
        2166           0 :       epsilon = 23.439_r8 - 4.e-7_r8*n
        2167           0 :       epsrad = epsilon*d2r
        2168             : !
        2169             : ! Right ascension:
        2170           0 :       alpha = atan2(cos(epsrad)*sinlam,cos(lamrad))*r2d
        2171             : !
        2172             : ! Declination:
        2173           0 :       delta = asin(sin(epsrad)*sinlam)*r2d
        2174             : !
        2175             : ! Subsolar latitude (output argument):
        2176           0 :       sbsllat = delta
        2177             : !
        2178             : ! Equation of time (degrees):
        2179           0 :       etdeg = l - alpha
        2180           0 :       nrot = nint(etdeg/360._r8)
        2181           0 :       etdeg = etdeg - dble(360*nrot)
        2182             : !
        2183             : ! Apparent time (degrees):
        2184             : ! Earth rotates one degree every 240 s.
        2185           0 :       aptime = ut/240._r8 + etdeg
        2186             : !
        2187             : ! Subsolar longitude (output argument):
        2188           0 :       sbsllon = 180._r8 - aptime
        2189           0 :       nrot = nint(sbsllon/360._r8)
        2190           0 :       sbsllon = sbsllon - dble(360*nrot)
        2191             :
        2192           0 : end subroutine apex_subsol
        2193             :
        2194             : !-----------------------------------------------------------------------
        2195           0 : subroutine solgmlon(xlat,xlon,colat,elon,mlon)
        2196             : !
        2197             : ! Compute geomagnetic longitude of the point with geocentric spherical
        2198             : !  latitude and longitude of XLAT and XLON, respectively.
        2199             : ! 940719 A. D. Richmond, NCAR
        2200             : !
        2201             : ! Algorithm:
        2202             : !   Use spherical coordinates.
        2203             : !   Let GP be geographic pole.
        2204             : !   Let GM be geomagnetic pole (colatitude COLAT, east longitude ELON).
        2205             : !   Let XLON be longitude of point P.
        2206             : !   Let TE be colatitude of point P.
        2207             : !   Let ANG be longitude angle from GM to P.
        2208             : !   Let TP be colatitude of GM.
        2209             : !   Let TF be arc length between GM and P.
        2210             : !   Let PA = MLON be geomagnetic longitude, i.e., Pi minus angle measured
        2211             : !     counterclockwise from arc GM-P to arc GM-GP.
        2212             : !   Then, using notation C=cos, S=sin, spherical-trigonometry formulas
        2213             : !     for the functions of the angles are as shown below.  Note: STFCPA,
        2214             : !     STFSPA are sin(TF) times cos(PA), sin(PA), respectively.
        2215             : !
        2216             : ! Input Args:
        2217             :   real(r8),intent(in)  :: xlat,xlon,colat,elon
        2218             : !
        2219             : ! Output Arg: Geomagnetic dipole longitude of the point (deg, -180. to 180.)
        2220             :   real(r8),intent(out) :: mlon
        2221             : !
        2222             : ! Local:
        2223             :   real(r8),parameter ::           &
        2224             :     rtod=5.72957795130823e1_r8,  &
        2225             :     dtor=1.745329251994330e-2_r8
        2226             :   real(r8) :: ctp,stp,ang,cang,sang,cte,ste,stfcpa,stfspa
        2227             :
        2228           0 :   ctp = cos(colat*dtor)
        2229           0 :   stp = sqrt(1._r8 - ctp*ctp)
        2230           0 :   ang = (xlon-elon)*dtor
        2231           0 :   cang = cos(ang)
        2232           0 :   sang = sin(ang)
        2233           0 :   cte = sin(xlat*dtor)
        2234           0 :   ste = sqrt(1._r8-cte*cte)
        2235           0 :   stfcpa = ste*ctp*cang - cte*stp
        2236           0 :   stfspa = sang*ste
        2237           0 :   mlon = atan2(stfspa,stfcpa)*rtod
        2238             :
        2239           0 : end subroutine solgmlon
        2240             : !-----------------------------------------------------------------------
        2241             : !-----------------------------------------------------------------------
        2242           0 : subroutine apex_magloctm (alon,sbsllat,sbsllon,clatp,polon,mlt)
        2243             :   !
        2244             :   !-----------------------------------------------------------------------
        2245             :   !  Computes magnetic local time from magnetic longitude, subsolar coordinates,
        2246             :   !   and geomagnetic pole coordinates.
        2247             :   !  950302 A. D. Richmond, NCAR
        2248             :   !  Algorithm:  MLT is calculated from the difference of the apex longitude,
        2249             :   !   alon, and the geomagnetic dipole longitude of the subsolar point.
        2250             :   !
        2251             :   !   Inputs:
        2252             :   !    alon    = apex magnetic longitude of the point (deg)
        2253             :   !    sbsllat = geographic latitude of subsolar point (degrees)
        2254             :   !    sbsllon = geographic longitude of subsolar point (degrees)
        2255             :   !    clatp   = Geocentric colatitude of geomagnetic dipole north pole (deg)
        2256             :   !    polon   = East longitude of geomagnetic dipole north pole (deg)
        2257             :   !
        2258             :   !   Output:
        2259             :   !    mlt (real) = magnetic local time for the apex longitude alon (hours)
        2260             :   !
        2261             :   !-----------------------------------------------------------------------
        2262             :   !
        2263             :   !------------------------------Arguments--------------------------------
        2264             :   !
        2265             :   REAL(r8) alon, sbsllat, sbsllon, clatp, polon, MLT
        2266             :   !
        2267             :   !---------------------------Local variables-----------------------------
        2268             :   !
        2269             :   real(r8) smlon
        2270             :   !
        2271             :   !-----------------------------------------------------------------------
        2272             :   !
        2273           0 :   call solgmlon (sbsllat,sbsllon,clatp,polon,smlon)
        2274           0 :   mlt = (alon - smlon)/15.0_r8 + 12.0_r8
        2275           0 :   if (mlt .ge. 24.0_r8) mlt = mlt - 24.0_r8
        2276           0 :   if (mlt .lt.   0._r8) mlt = mlt + 24.0_r8
        2277           0 :   return
        2278             : end subroutine apex_magloctm
        2279             :
        2280             : !-----------------------------------------------------------------------
        2281             : end module apex

Generated by: LCOV version 1.14