From 3cbf5576da9c9a443d077c4b4f4e2441fb192509 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sat, 10 Feb 2024 13:11:10 -0600 Subject: [PATCH] restored shellc as a separate subroutine --- src/shellig.f90 | 65 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 54 insertions(+), 11 deletions(-) diff --git a/src/shellig.f90 b/src/shellig.f90 index cc8e690..86a6a3d 100644 --- a/src/shellig.f90 +++ b/src/shellig.f90 @@ -247,6 +247,34 @@ subroutine findb0(me,stps,bdel,value,bequ,rr0) end subroutine findb0 +!***************************************************************************************** +!> +! Wrapper to [[shellg]] to be used with cartesian coordinates. +! +!@note In the original code, this was an ENTRY point in [[shellg]] and didn't +! include all the outputs. + + subroutine shellc(me,v,dimo,fl,icode,b0) + + class(shellig_type),intent(inout) :: me + real(wp),dimension(3),intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km) + !! * x-axis pointing to equator at 0 longitude + !! * y-axis pointing to equator at 90 long. + !! * z-axis pointing to north pole + real(wp),intent(in) :: dimo !! dipol moment in gauss (normalized to earth radius) + real(wp),intent(out) :: fl !! l-value + integer,intent(out) :: icode !! * =1 normal completion + !! * =2 unphysical conjugate point (fl meaningless) + !! * =3 shell parameter greater than limit up to + !! which accurate calculation is required; + !! approximation is used. + real(wp),intent(out) :: b0 !! magnetic field strength in gauss + real(wp) :: glat,glon,alt !! not used + + call me%shellg(glat,glon,alt,dimo,fl,icode,b0,v) + + end subroutine shellc + !***************************************************************************************** !> ! calculates l-value for specified geodaetic coordinates, altitude @@ -262,7 +290,7 @@ end subroutine findb0 ! - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/ ! - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990 - subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0) + subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0,v) class(shellig_type),intent(inout) :: me real(wp),intent(in) :: glat !! geodetic latitude in degrees (north) @@ -276,12 +304,20 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0) !! which accurate calculation is required; !! approximation is used. real(wp),intent(out) :: b0 !! magnetic field strength in gauss + real(wp),dimension(3),intent(in),optional :: v !! cartesian coordinates in earth radii (6371.2 km) + !! + !! * x-axis pointing to equator at 0 longitude + !! * y-axis pointing to equator at 90 long. + !! * z-axis pointing to north pole + !! + !! If this argument is present, it is used + !! instead of glat,glon,alt. See [[shellc]]. real(wp) :: arg1 , arg2 , bequ , bq1 , bq2 , bq3 , c0 , c1 , c2 , c3 , & ct , d , d0 , d1 , d2, dimob0 , e0 , e1 , e2 , ff , fi , gg , & hli , oradik , oterm , p(8,100) , r , r1 , r2 , r3 , r3h , radik , & rlat , rlon , rq , st , step12 , step2 , & - stp , t , term , v(3) , xx , z , zq , zz + stp , t , term , xx , z , zq , zz integer :: i , iequ , n real(wp),parameter :: rmin = 0.05_wp !! boundaries for identification of `icode=2 and 3` @@ -289,15 +325,22 @@ subroutine shellg(me,glat,glon,alt,dimo,fl,icode,b0) integer,parameter :: max_loop_index = 100 ! 3333 <--- jw : original code had 3333 ... was this a bug ???? bequ = 1.0e10_wp - rlat = glat*umr - ct = sin(rlat) - st = cos(rlat) - d = sqrt(aquad-(aquad-bquad)*ct*ct) - me%xi(1) = (alt+aquad/d)*st/era - me%xi(3) = (alt+bquad/d)*ct/era - rlon = glon*umr - me%xi(2) = me%xi(1)*sin(rlon) - me%xi(1) = me%xi(1)*cos(rlon) + + if (present(v)) then + me%xi(1) = v(1) + me%xi(2) = v(2) + me%xi(3) = v(3) + else + rlat = glat*umr + rlon = glon*umr + ct = sin(rlat) + st = cos(rlat) + d = sqrt(aquad-(aquad-bquad)*ct*ct) + me%xi(1) = (alt+aquad/d)*st/era + me%xi(3) = (alt+bquad/d)*ct/era + me%xi(2) = me%xi(1)*sin(rlon) + me%xi(1) = me%xi(1)*cos(rlon) + end if !*****convert to dipol-oriented co-ordinates rq = 1.0_wp/(me%xi(1)*me%xi(1)+me%xi(2)*me%xi(2)+me%xi(3)*me%xi(3))