subroutine cread (in,out,ty,l,maxln,a) ! reads schmidt coeffs (with true second derivatives) ! and sets up for gfield and/or sphrc ! mod 5/15/83 to make less work of time update (j. cain) ! only updates after l=0 to maximum nonzero gt or ht value ! expanded coefficient format: 2i3,2f16.8,f,4f10.5 ! also reads from unit 2, not 12 parameter (md=101) common /coeff/tg(md,md),const(md,md),fm(md),fn(md),ge(8) dimension g(md,md),gt(md,md),gtt(md,md) character (len=72) :: ai integer out if (l /= 0) then minn = 2 call con maxn = 0 maxnt = 0 read (in,9) ai read (in,*) a,tzero,tf,tl do read (in,10) ln,lm,gnm,hnm,gtnm,htnm,gttnm,httnm n = ln + 1 m = lm + 1 if (lm*ln == 0) minn = 1 if (ln < 0) exit maxn = max(n,maxn) if (gtnm /= 0. .or. htnm /= 0.) maxnt = max(n,maxnt) g(n,m) = gnm gt(n,m) = gtnm gtt(n,m) = gttnm if (lm == 0) cycle g(lm,n) = hnm gt(lm,n) = htnm gtt(lm,n) = httnm enddo maxln = maxn - 1 read (in,16) ge write (out,11) ai,maxln,maxnt-1,a,tzero,tf,tl if (l <= 1) then write (out,12) do n = minn, maxn ln = n - 1 do m = 1, n lm = m - 1 if (m == 1) then write (out,14) ln,lm,g(n,m),gt(n,m),gtt(n,m) else write (out,13) ln,lm,g(n,m),g(lm,n),gt(n,m),gt(lm,n),gtt(n,m),gtt(lm,n) endif enddo enddo write (out,44) ge endif do n = 1, maxn do m = 1, maxn call convrt(g(n,m),n,m,1) tg(n,m) = g(n,m) call convrt(gt(n,m),n,m,1) gtt(n,m) = 0.5*gtt(n,m) call convrt(gtt(n,m),n,m,1) enddo enddo if(ty > tl .or. ty < tf) write (out,15) ty l = 0 endif t = ty - tzero do n = 1, maxnt do m = 1, maxnt tg(n,m) = g(n,m) + t*(gt(n,m)+t*gtt(n,m)) enddo enddo return ! two formats for coeffs ! 10 format (2i3,6f11.5) 10 format (2i3,2f16.8,4f10.5) 9 format (a) 11 format ('0',a,/2x,'mxln=',i3,2x,'mxlnt=',i3,2x,'a=',f6.1,2x,'tzero=',f7.1,2x,'tf=',f7.1,2x,'tl=',f7.1) 12 format ('0 n m',8x,'g',10x,'h',10x,'gt',9x,'ht',8x,'gtt',8x,'htt') 13 format (2i3,2f11.2,2f11.3,2f11.4) 14 format (2i3,f11.2,11x,f11.3,11x,f11.4) 15 format ('0**warning,',f7.1,' is outside time limits') 16 format (6x,8f11.4) 44 format (' first externals',/6x,8f11.4) end subroutine convrt (g,i,l,k) ! k=1 converts schmidt to gauss,k>1 converts gauss to schmidt parameter (md=101) dimension s(md,md) logical next save s data next/.false./ if (.not. next) then next = .true. s(1,1) = -1. do n = 2, md s(n,1) = s(n-1,1)*real(2*n-3)/real(n-1) s(1,n) = 0. j = 2 do m = 2, n s(n,m) = s(n,m-1)*sqrt(real((n-m+1)*j)/real(n+m-2)) s(m-1,n) = s(n,m) j = 1 enddo enddo endif if (k > 1) then g = g/s(i,l) else g = g*s(i,l) endif return end subroutine con ! sets up constants for use by sphrc (normally called by one ! of the subroutines that sets up the coefficients tg) parameter (md=101) common/coeff/tg(md,md),c(md,md),fm(md),fn(md),ge(8) ! data c/3721*0./ ! assume computer core preset to zero or needed fm(1) = 0. do n = 2, md fm(n) = n-1 fn(n) = n do m = 1, n c(n,m) = real( (n-2)**2 - (m-1)**2 )/real( (2*n-3)*(2*n-5) ) enddo enddo return end subroutine sphrc ! assume computer core is initiallized to zero ! otherwise should zero p,dp,sp arrays ! this version best if frequent dipole only computation. ! otherwise would be better to put n=2 in main loop parameter (md=101) common /gcom/ st,ct,sph,cph,aor,bt,bp,br,nmax,sind,cosd common /coeff/ g(md,md),const(md,md),fm(md),fn(md),ge(8) dimension p(md,md),dp(md,md),sp(md),cp(md) dimension pp(0:md,0:md),dpp(0:md,0:md),s(0:md,0:md),q(0:md,0:md) ar2 = aor**2 ar = ar2*aor gc = g(2,2)*cph + g(1,2)*sph br = -ar2*g(1,1) - 2.*ar*(g(2,1)*ct+gc*st) bt = ar*(gc*ct-g(2,1)*st) bp = ar*(g(1,2)*cph-g(2,2)*sph) if (nmax /= 2) then bp = bp*st p(2,1) = ct dp(2,1) = -st p(2,2) = st dp(2,2) = ct sp(2) = sph cp(2) = cph p(1,1) = 1. do n = 3, nmax sp(n) = sph*cp(n-1) + cph*sp(n-1) cp(n) = cph*cp(n-1) - sph*sp(n-1) dp(n,1) = ct*dp(n-1,1) - st*p(n-1,1) - const(n,1)*dp(n-2,1) stm = g(n,1)*dp(n,1) spm = 0. p(n,1) = ct*p(n-1,1) - const(n,1)*p(n-2,1) srm = -g(n,1)*p(n,1) p(n,n) = st*p(n-1,n-1) dp(n,n) = fm(n)*ct*p(n-1,n-1) do m = 2, n if (n /= m) then p(n,m) = ct*p(n-1,m) - const(n,m)*p(n-2,m) dp(n,m) = ct*dp(n-1,m) - st*p(n-1,m) - const(n,m)*dp(n-2,m) endif gc = g(n,m)*cp(m) + g(m-1,n)*sp(m) stm = stm + gc*dp(n,m) spm = spm + (g(m-1,n)*cp(m)-g(n,m)*sp(m))*fm(m)*p(n,m) srm = srm - gc*p(n,m) enddo ar = aor*ar bt = bt + stm*ar bp = bp + spm*ar br = br + srm*fn(n)*ar enddo if(st == 0.) st = 1.0e-9 bp = bp/st endif ! assume first order external usually non-zero roa = 1./aor c2ph = 2.*cph*cph - 1. s2ph = 2.*sph*cph q(1,0) = ge(1) q(1,1) = ge(2) s(1,1) = ge(3) q(2,0) = ge(4) q(2,1) = ge(5) s(2,1) = ge(6) q(2,2) = ge(7) s(2,2) = ge(8) pp(1,0) = p(2,1) pp(1,1) = p(2,2) pp(2,0) = p(3,1) pp(2,1) = p(3,2) pp(2,2) = p(3,3) dpp(1,0) = dp(2,1) dpp(1,1) = dp(2,2) dpp(2,0) = dp(3,1) dpp(2,1) = dp(3,2) dpp(2,2) = dp(3,3) e_r = q(1,0)*pp(1,0)+(q(1,1)*cph+s(1,1)*sph)*pp(1,1) & + (q(2,0)*pp(2,0)+(q(2,1)*cph+s(2,1)*sph)*pp(2,1)+(q(2,2)*c2ph+s(2,2)*s2ph)*pp(2,2))*2.0*roa e_t = q(1,0)*dpp(1,0)+(q(1,1)*cph+s(1,1)*sph)*dpp(1,1) & + (q(2,0)*dpp(2,0)+(q(2,1)*cph+s(2,1)*sph)*dpp(2,1)+(q(2,2)*c2ph+s(2,2)*s2ph)*dpp(2,2))*roa e_p = (s(1,1)*cph-q(1,1)*sph) & + ((s(2,1)*cph-q(2,1)*sph)*pp(2,1)+(s(2,2)*c2ph-q(2,2)*s2ph)*2.0*pp(2,2))*roa/st br = br - e_r bt = bt - e_t bp = bp - e_p return end subroutine gfield(dlat,dlong,alt,mln,a,x,y,z,f) common/gcom/st,ct,sph,cph,aor,bt,bp,br,nmax,sind,cosd save re,a2,a4,a2b2,b2,a4b4 data re,tlast,rp/2*0.,3374.9/ if (re /= 3396.9) then re = 3396.9 a2 = re**2 a4 = a2**2 ! flat = 1.-1./298.25 flat=rp/re b2 = (re*flat)**2 a2b2 = a2*(1.-flat**2) a4b4 = a4*(1.-flat**4) drad = atan(1.)/45. end if sinla = sin(dlat*drad) rlong = dlong*drad cph = cos(rlong) sph = sin(rlong) sinla2 = sinla**2 cosla2 = 1. - sinla**2 den2 = a2 - a2b2*sinla**2 den = sqrt(den2) fac = ((alt*den+a2)/(alt*den+b2))**2 r = sqrt(alt*(alt+2.*den)+(a4-a4b4*sinla2)/den2) ct = sinla/sqrt(fac*cosla2+sinla2) st = sqrt(1.-ct**2) nmax = mln + 1 aor = a/r call sphrc f = sqrt(bt**2+bp**2+br**2) ! transforms to geodetic directions sind = sinla*st-sqrt(cosla2)*ct cosd = sqrt(1.-sind**2) x = -bt*cosd-br*sind y = bp z = bt*sind-br*cosd return end