Skip to content

Commit

Permalink
eliminated the entry routines from shellig
Browse files Browse the repository at this point in the history 
added a new igrfc routine (untested)
  • Loading branch information
@jacobwilliams
jacobwilliams committed Feb 10, 2024
1 parent 3cbf557 commit 39e720d
Showing 1 changed file with 233 additions and 102 deletions.
335 changes: 233 additions & 102 deletions src/shellig.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,13 +69,13 @@ module shellig_module
contains
private

procedure,public :: igrf
procedure,public :: igrf, igrfc

procedure, public :: feldcof
procedure, public :: feldg
procedure, public :: shellg
procedure, public :: feldg, feldc
procedure, public :: shellg, shellc
procedure, public :: findb0
procedure :: stoer, getshc, intershc, extrashc
procedure :: stoer, getshc, intershc, extrashc, feldi
procedure,public :: set_data_file_dir, get_data_file_dir

end type shellig_type
Expand Down Expand Up @@ -151,6 +151,51 @@ subroutine igrf(me,lon,lat,height,year,xl,bbx)
end subroutine igrf
!*****************************************************************************************

!*****************************************************************************************
!>
! Alternate version of [[igrf]] for cartesian coordinates.

subroutine igrfc(me,v,year,xl,bbx)

class(shellig_type),intent(inout) :: me
real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km)
!! x-axis pointing to equator at 0 longitude
!! y-axis pointing to equator at 90 long.
!! z-axis pointing to north pole
real(wp),intent(in) :: year !! decimal year for which geomagnetic field is to
!! be calculated (e.g.:1995.5 for day 185 of 1995)
real(wp),intent(out) :: xl !! l-value
real(wp),intent(out) :: bbx !! b_total / b_equatorial ratio

real(wp) :: bab1 , babs , bdel , beq , bequ , dimo , rr0
integer :: icode
logical :: val
real(wp),dimension(3) :: b

real(wp),parameter :: stps = 0.05_wp

! JW : do we need to reset some or all of these ?
me%sp = 0.0_wp
me%xi = 0.0_wp
me%h = 0.0_wp
me%step = 0.20_wp
me%steq = 0.03_wp

call me%feldcof(year,dimo)
call me%feldc(v,b)
call me%shellc(v,dimo,xl,icode,bab1)

bequ = dimo/(xl*xl*xl)
if ( icode==1 ) then
bdel = 1.0e-3_wp
call me%findb0(stps,bdel,val,beq,rr0)
if ( val ) bequ = beq
endif
bbx = babs/bequ

end subroutine igrfc
!*****************************************************************************************

!*****************************************************************************************
!>
subroutine findb0(me,stps,bdel,value,bequ,rr0)
Expand Down Expand Up @@ -516,7 +561,7 @@ end subroutine shellg

!*****************************************************************************************
!>
! subroutine used for field line tracing in [[shellg]]
! subroutine used for field line tracing in [[shellg]].
! calls entry point [[feldi]] in geomagnetic field subroutine [[feldg]]

subroutine stoer(me,p,bq,r)
Expand Down Expand Up @@ -545,7 +590,7 @@ subroutine stoer(me,p,bq,r)
! Changed from CALL FELDI(XI,H); XI, H are in COMMON block; results
! are the same; dkb Feb 1998.
! JW : feb 2024 : xi, h now class variables.
CALL feldi(me)
call me%feldi()
q = me%H(1)/rq
dx = me%H(3) + me%H(3) + q*me%Xi(1)
dy = me%H(4) + me%H(4) + q*me%Xi(2)
Expand All @@ -562,7 +607,7 @@ subroutine stoer(me,p,bq,r)
Bq = dsq*rq*rq
P(6) = sqrt(dsq/(rq+3.0_wp*zm*zm))
P(7) = P(6)*(rq+zm*zm)/(rq*dzm)
END subroutine stoer
end subroutine stoer

!*****************************************************************************************
!>
Expand All @@ -576,8 +621,11 @@ END subroutine stoer
! * changes (d. bilitza, nov 87):
! - field coefficients in binary data files instead of block data
! - calculates dipol moment
!
!@note In the original code, [[feldc] and [[feldi]] were
! ENTRY points to this routine

subroutine feldg(me,glat,glon,alt,bnorth,beast,bdown,babs)
subroutine feldg(me,glat,glon,alt,bnorth,beast,bdown,babs)

class(shellig_type),intent(inout) :: me
real(wp),intent(in) :: glat !! geodetic latitude in degrees (north)
Expand All @@ -589,104 +637,187 @@ subroutine feldg(me,glat,glon,alt,bnorth,beast,bdown,babs)
!! and downward.
real(wp),intent(out) :: Babs !! magnetic field strength in gauss

real(wp) :: b(3) , brho , bxxx , &
byyy , bzzz , cp , ct , d , f , rho , &
rlat , rlon , rq , s , sp , st , t , v(3) , x , xxx , &
y , yyy , z , zzz
integer :: i , ih , ihmax , il , imax , is , k , last , m

!-- is records entry point
!
!*****entry point feldg to be used with geodetic co-ordinates
is=1
rlat=glat*umr
ct=sin(rlat)
st=cos(rlat)
d=sqrt(aquad-(aquad-bquad)*ct*ct)
rlon=glon*umr
cp=cos(rlon)
sp=sin(rlon)
zzz=(alt+bquad/d)*ct/era
rho=(alt+aquad/d)*st/era
xxx=rho*cp
yyy=rho*sp
goto 10

!*****entry point feldc to be used with cartesian co-ordinates
! v(3) cartesian coordinates in earth radii (6371.2 km)
! x-axis pointing to equator at 0 longitude
! y-axis pointing to equator at 90 long.
! z-axis pointing to north pole
entry feldc(me,v,b)
is=2
xxx=v(1)
yyy=v(2)
zzz=v(3)

10 rq=1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz)
me%xi(1)=xxx*rq
me%xi(2)=yyy*rq
me%xi(3)=zzz*rq
goto 20

!*****entry point feldi used for l computation
entry feldi(me)
is=3
20 ihmax=me%nmax*me%nmax+1
last=ihmax+me%nmax+me%nmax
imax=me%nmax+me%nmax-1
do i=ihmax,last
me%h(i)=me%g(i)
end do
do k=1,3,2
i=imax
ih=ihmax
do
il=ih-i
f=2.0_wp/real(i-k+2, wp)
x=me%xi(1)*f
y=me%xi(2)*f
z=me%xi(3)*(f+f)
i=i-2
if ((i-1)>=0) then
if ((i-1)>0) then
do m=3,i,2
me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-&
me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2))
me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-&
me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1))
end do
end if
me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih))
me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih))
end if
me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2))
ih=il
if (i<k) exit
end do
end do

if (is==3) return

s=0.5_wp*me%h(1)+2.0_wp*(me%h(2)*me%xi(3)+me%h(3)*me%xi(1)+me%h(4)*me%xi(2))
t=(rq+rq)*sqrt(rq)
bxxx=t*(me%h(3)-s*xxx)
byyy=t*(me%h(4)-s*yyy)
bzzz=t*(me%h(2)-s*zzz)
if (is==2) then
b(1)=bxxx
b(2)=byyy
b(3)=bzzz
else
babs=sqrt(bxxx*bxxx+byyy*byyy+bzzz*bzzz)
beast=byyy*cp-bxxx*sp
brho=byyy*sp+bxxx*cp
bnorth=bzzz*st-brho*ct
bdown=-bzzz*ct-brho*st
real(wp) :: brho , bxxx , byyy , bzzz , cp , ct , d , f , rho , &
rlat , rlon , rq , s , sp , st , t , &
x , xxx , y , yyy , z , zzz
integer :: i , ih , ihmax , il , imax , k , last , m

rlat=glat*umr
ct=sin(rlat)
st=cos(rlat)
d=sqrt(aquad-(aquad-bquad)*ct*ct)
rlon=glon*umr
cp=cos(rlon)
sp=sin(rlon)
zzz=(alt+bquad/d)*ct/era
rho=(alt+aquad/d)*st/era
xxx=rho*cp
yyy=rho*sp

rq=1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz)
me%xi = [xxx,yyy,zzz] * rq

ihmax=me%nmax*me%nmax+1
last=ihmax+me%nmax+me%nmax
imax=me%nmax+me%nmax-1
do i=ihmax,last
me%h(i)=me%g(i)
end do
do k=1,3,2
i=imax
ih=ihmax
do
il=ih-i
f=2.0_wp/real(i-k+2, wp)
x=me%xi(1)*f
y=me%xi(2)*f
z=me%xi(3)*(f+f)
i=i-2
if ((i-1)>=0) then
if ((i-1)>0) then
do m=3,i,2
me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-&
me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2))
me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-&
me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1))
end do
end if
me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih))
me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih))
end if
me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2))
ih=il
if (i<k) exit
end do
end do

s=0.5_wp*me%h(1)+2.0_wp*(me%h(2)*me%xi(3)+me%h(3)*me%xi(1)+me%h(4)*me%xi(2))
t=(rq+rq)*sqrt(rq)
bxxx=t*(me%h(3)-s*xxx)
byyy=t*(me%h(4)-s*yyy)
bzzz=t*(me%h(2)-s*zzz)

babs=sqrt(bxxx*bxxx+byyy*byyy+bzzz*bzzz)
beast=byyy*cp-bxxx*sp
brho=byyy*sp+bxxx*cp
bnorth=bzzz*st-brho*ct
bdown=-bzzz*ct-brho*st

end subroutine feldg

!*****************************************************************************************
!>
! Alternate version of [[feldg]] to be used with cartesian coordinates

subroutine feldc(me,v,b)

class(shellig_type),intent(inout) :: me
real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km)
!! x-axis pointing to equator at 0 longitude
!! y-axis pointing to equator at 90 long.
!! z-axis pointing to north pole
real(wp),intent(out) :: b(3) !! field components

real(wp) :: f , rq , s , t , x , xxx , y , yyy , z , zzz
integer :: i , ih , ihmax , il , imax , k , last , m

xxx=v(1)
yyy=v(2)
zzz=v(3)

rq=1.0_wp/(xxx*xxx+yyy*yyy+zzz*zzz)
me%xi = [xxx,yyy,zzz] * rq

ihmax=me%nmax*me%nmax+1
last=ihmax+me%nmax+me%nmax
imax=me%nmax+me%nmax-1
do i=ihmax,last
me%h(i)=me%g(i)
end do
do k=1,3,2
i=imax
ih=ihmax
do
il=ih-i
f=2.0_wp/real(i-k+2, wp)
x=me%xi(1)*f
y=me%xi(2)*f
z=me%xi(3)*(f+f)
i=i-2
if ((i-1)>=0) then
if ((i-1)>0) then
do m=3,i,2
me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-&
me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2))
me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-&
me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1))
end do
end if
me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih))
me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih))
end if
me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2))
ih=il
if (i<k) exit
end do
end do

s=0.5_wp*me%h(1)+2.0_wp*(me%h(2)*me%xi(3)+me%h(3)*me%xi(1)+me%h(4)*me%xi(2))
t=(rq+rq)*sqrt(rq)

b(1)=t*(me%h(3)-s*xxx)
b(2)=t*(me%h(4)-s*yyy)
b(3)=t*(me%h(2)-s*zzz)

end subroutine feldc

!*****************************************************************************************
!>
! Used for `l` computation.

subroutine feldi(me)

class(shellig_type),intent(inout) :: me

real(wp) :: f , x , y , z
integer :: i , ih , ihmax , il , imax , k , last , m

ihmax=me%nmax*me%nmax+1
last=ihmax+me%nmax+me%nmax
imax=me%nmax+me%nmax-1
do i=ihmax,last
me%h(i)=me%g(i)
end do
do k=1,3,2
i=imax
ih=ihmax
do
il=ih-i
f=2.0_wp/real(i-k+2, wp)
x=me%xi(1)*f
y=me%xi(2)*f
z=me%xi(3)*(f+f)
i=i-2
if ((i-1)>=0) then
if ((i-1)>0) then
do m=3,i,2
me%h(il+m+1)=me%g(il+m+1)+z*me%h(ih+m+1)+x*(me%h(ih+m+3)-&
me%h(ih+m-1))-y*(me%h(ih+m+2)+me%h(ih+m-2))
me%h(il+m)=me%g(il+m)+z*me%h(ih+m)+x*(me%h(ih+m+2)-&
me%h(ih+m-2))+y*(me%h(ih+m+3)+me%h(ih+m-1))
end do
end if
me%h(il+2)=me%g(il+2)+z*me%h(ih+2)+x*me%h(ih+4)-y*(me%h(ih+3)+me%h(ih))
me%h(il+1)=me%g(il+1)+z*me%h(ih+1)+y*me%h(ih+4)+x*(me%h(ih+3)-me%h(ih))
end if
me%h(il)=me%g(il)+z*me%h(ih)+2.0_wp*(x*me%h(ih+1)+y*me%h(ih+2))
ih=il
if (i<k) exit
end do
end do

end subroutine feldi

!*****************************************************************************************
!>
! Determines coefficients and dipol moment from IGRF models
Expand Down

0 comments on commit 39e720d

Please sign in to comment.