restored shellc as a separate subroutine

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,28 +304,43 @@ 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`
    real(wp),parameter :: rmax = 1.01_wp !! boundaries for identification of `icode=2 and 3`
    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))