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 16128 : 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 16128 : if (igrf_set) return
231 :
232 : !----------------------------------------------------------------------
233 : ! ... open the netcdf file
234 : !----------------------------------------------------------------------
235 1536 : call getfil(coefs_file, locfn, 0)
236 1536 : call cam_pio_openfile( ncid, trim(locfn), PIO_NOWRITE)
237 :
238 : !----------------------------------------------------------------------
239 : ! ... read the snoe dimensions
240 : !----------------------------------------------------------------------
241 1536 : ierr = pio_inq_dimid( ncid, 'n1', dim_id )
242 1536 : ierr = pio_inq_dimlen( ncid, dim_id, n1 )
243 1536 : ierr = pio_inq_dimid( ncid, 'ncn1', dim_id )
244 1536 : ierr = pio_inq_dimlen( ncid, dim_id, ncn1 )
245 :
246 1536 : ierr = pio_inq_dimid( ncid, 'n2', dim_id )
247 1536 : ierr = pio_inq_dimlen( ncid, dim_id, n2 )
248 1536 : ierr = pio_inq_dimid( ncid, 'ncn2', dim_id )
249 1536 : ierr = pio_inq_dimlen( ncid, dim_id, ncn2 )
250 :
251 10752 : allocate( g1(n1,ncn1), g2(n2,ncn2) )
252 :
253 1536 : ierr = pio_inq_varid( ncid, 'g1', var_id )
254 1536 : ierr = pio_get_var( ncid, var_id, g1 )
255 :
256 1536 : ierr = pio_inq_varid( ncid, 'g2', var_id )
257 1536 : ierr = pio_get_var( ncid, var_id, g2 )
258 :
259 1536 : ierr = pio_inq_varid( ncid, 'era1_year', var_id )
260 1536 : ierr = pio_get_var( ncid, var_id, year1 )
261 :
262 1536 : ierr = pio_inq_varid( ncid, 'era2_year', var_id )
263 1536 : ierr = pio_get_var( ncid, var_id, year2 )
264 :
265 1536 : call pio_closefile(ncid)
266 :
267 1536 : apex_beg_yr = year1
268 1536 : apex_end_yr = year2+5*(ncn2-1)
269 :
270 1536 : igrf_set = .true.
271 :
272 17664 : end subroutine apex_set_igrf
273 :
274 : !-----------------------------------------------------------------------
275 1536 : 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 1536 : 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 1536 : rtd = 45._r8/atan(1._r8)
308 1536 : 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 1536 : 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 1536 : if (.not.allocated(xarray)) then
323 7680 : allocate(xarray(nlat,nlon,nalt),stat=istat)
324 1536 : if (istat /= 0) call endrun( 'allocate xarray' )
325 201841152 : xarray = 0._r8
326 : endif
327 1536 : if (.not.allocated(yarray)) then
328 7680 : allocate(yarray(nlat,nlon,nalt),stat=istat)
329 1536 : if (istat /= 0) call endrun( 'allocate yarray' )
330 201841152 : yarray = 0._r8
331 : endif
332 1536 : if (.not.allocated(zarray)) then
333 7680 : allocate(zarray(nlat,nlon,nalt),stat=istat)
334 1536 : if (istat /= 0) call endrun( 'allocate zarray' )
335 201841152 : zarray = 0._r8
336 : endif
337 1536 : if (.not.allocated(varray)) then
338 7680 : allocate(varray(nlat,nlon,nalt),stat=istat)
339 1536 : if (istat /= 0) call endrun( 'allocate varray' )
340 201841152 : 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 1536 : nglon=nlon ; nglat=nlat ; ngalt=nalt
347 4608 : if (.not.allocated(geolat)) allocate(geolat(nglat),stat=istat)
348 4608 : if (.not.allocated(geolon)) allocate(geolon(nglon),stat=istat)
349 4608 : if (.not.allocated(geoalt)) allocate(geoalt(ngalt),stat=istat)
350 279552 : geolat(:) = gplat(:)
351 556032 : geolon(:) = gplon(:)
352 4608 : geoalt(:) = gpalt(:)
353 : !
354 : ! Set coefficients gb,gv (module data) for requested year:
355 : !
356 1536 : 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 1536 : call apex_dypol(colat,elon,vp)
362 :
363 1536 : ctp = cos(colat*dtr)
364 1536 : stp = sin(colat*dtr)
365 :
366 1536 : reqore = req/re
367 1536 : rqorm1 = reqore-1._r8
368 :
369 279552 : do j=1,nlat
370 278016 : ct = sin(gplat(j)*dtr)
371 278016 : st = cos(gplat(j)*dtr)
372 278016 : kpol = 0
373 278016 : if (abs(gplat(j)) > pola) kpol = 1
374 100643328 : do i=1,nlon
375 100363776 : if (kpol==1.and.i > 1) then
376 3317760 : xarray(j,i,:) = xarray(j,1,:)
377 3317760 : yarray(j,i,:) = yarray(j,1,:)
378 3317760 : zarray(j,i,:) = zarray(j,1,:)
379 3317760 : varray(j,i,:) = varray(j,1,:)
380 : cycle
381 : endif
382 99257856 : cp = cos((gplon(i)-elon)*dtr)
383 99257856 : 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 99257856 : ctm = ctp*ct + stp*st*cp
391 99257856 : stmcpm = st*ctp*cp - ct*stp
392 99257856 : stmspm = st*sp
393 298051584 : do k=1,nalt
394 198515712 : call apex_sub(date,gplat(j),gplon(i),gpalt(k),&
395 198515712 : aht,alat,phia,bmag,xmag,ymag,zdown,vmp)
396 :
397 198515712 : vnor = vmp/vp
398 198515712 : rp = 1._r8 + gpalt(k)/re
399 198515712 : varray(j,i,k) = vnor*rp*rp + ctm
400 198515712 : reqam1 = req*(aht-1._r8)
401 198515712 : slp = sqrt(max(reqam1-gpalt(k),0._r8)/(reqam1+re))
402 : !
403 : ! Reverse sign of slp in southern magnetic hemisphere
404 : !
405 198515712 : if (zdown.lt.0._r8) slp = -slp
406 198515712 : clp = sqrt (rp/(reqore*aht-rqorm1))
407 198515712 : phiar = phia*dtr
408 198515712 : xarray(j,i,k) = clp*cos (phiar) - stmcpm
409 198515712 : yarray(j,i,k) = clp*sin (phiar) - stmspm
410 496289280 : 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 1536 : glatmx = max( glatlim,gplat(nlat-2))
420 1536 : glatmn = min(-glatlim,gplat(2))
421 :
422 16128 : end subroutine apex_mka
423 : !-----------------------------------------------------------------------
424 336384 : 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 112128 : ier = 0
475 112128 : 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 112128 : dfxdh ,dfydh ,dfzdh ,dfvdh, ier)
482 :
483 112128 : 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 112128 : 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 112128 : 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 112128 : if (glat > glatmx .or. glat < glatmn) then
503 1152 : glatx = glatmx
504 1152 : 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 1152 : dfvdln,dmxdh,dmydh,dmzdh,dmvdh, ier)
510 :
511 : call adpl(glatx,glonloc,cth,sth,fxdum,fydum,fzdum,fvdum, &
512 1152 : dmxdth,dmydth,dmzdth,dmvdth,dfxdln,dfydln,dfzdln,dfvdln)
513 :
514 : call grapxyzv(alt,cth,sth,dfxdln,dfydln,dfzdln,dfvdln, &
515 1152 : gradx,grady,gradz,gradv)
516 : endif
517 :
518 : call gradlpv(hr,alt,fx,fy,fz,fv,gradx,grady,gradz,gradv, &
519 112128 : 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 112128 : bmag,sim,si,f,d,w,bhat,d1,d2,d3,e1,e2,e3,f1,f2)
523 :
524 112128 : be3 = bmag/d
525 112128 : 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 112128 : 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 112128 : d2 = 40680925.e0_r8 - 272340.e0_r8*cth*cth
728 112128 : d = sqrt(d2)
729 112128 : rho = sth*(alt + 40680925.e0_r8/d)
730 112128 : dddthod = 272340.e0_r8*cth*sth/d2
731 112128 : drhodth = alt*cth + (40680925.e0_r8/d)*(cth-sth*dddthod)
732 112128 : dzetdth =-alt*sth - (40408585.e0_r8/d)*(sth+cth*dddthod)
733 112128 : ddisdth = sqrt(drhodth*drhodth + dzetdth*dzetdth)
734 :
735 112128 : gradx(1) = dfxdln/rho
736 112128 : grady(1) = dfydln/rho
737 112128 : gradz(1) = dfzdln/rho
738 112128 : gradv(1) = dfvdln/rho
739 :
740 112128 : gradx(2) = -dfxdth/ddisdth
741 112128 : grady(2) = -dfydth/ddisdth
742 112128 : gradz(2) = -dfzdth/ddisdth
743 112128 : gradv(2) = -dfvdth/ddisdth
744 :
745 112128 : gradx(3) = dfxdh
746 112128 : grady(3) = dfydh
747 112128 : gradz(3) = dfzdh
748 112128 : gradv(3) = dfvdh
749 :
750 112128 : end subroutine gradxyzv
751 : !-----------------------------------------------------------------------
752 1152 : 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 1152 : d2 = 40680925.e0_r8 - 272340.e0_r8*cth*cth
770 1152 : d = sqrt(d2)
771 1152 : rho = sth*(alt + 40680925.e0_r8/d)
772 :
773 1152 : gradx(1) = dfxdln/rho
774 1152 : grady(1) = dfydln/rho
775 1152 : gradz(1) = dfzdln/rho
776 1152 : gradv(1) = dfvdln/rho
777 :
778 1152 : end subroutine grapxyzv
779 : !-----------------------------------------------------------------------
780 112128 : 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 112128 : xlatm=0._r8 ; xlonm=0._r8 ; vmp=0._r8 ; grclm=0._r8 ; clmgrp=0._r8
815 112128 : rgrlp = 0._r8 ; b=0._r8 ; clm=0._r8 ; r3_2=0._r8 ; qdlat=0._r8
816 :
817 112128 : rr = re + hr
818 112128 : r = re + alt
819 112128 : rn = r/re
820 112128 : sqrror = sqrt(rr/r)
821 112128 : r3_2 = 1._r8/sqrror/sqrror/sqrror
822 112128 : xlonm = atan2(fy,fx)
823 112128 : cpm = cos(xlonm)
824 112128 : spm = sin(xlonm)
825 112128 : xlonm = rtd*xlonm ! output
826 112128 : bo = vp*1.e6_r8 ! vp is module data; 1.e6 converts T to nT and km-1 to m-1
827 112128 : rn2 = rn*rn
828 112128 : vmp = vp*fv/rn2 ! output
829 112128 : b(1) = -bo*gradv(1)/rn2
830 112128 : b(2) = -bo*gradv(2)/rn2
831 112128 : b(3) = -bo*(gradv(3)-2._r8*fv/r)/rn2
832 :
833 112128 : x2py2 = fx*fx + fy*fy
834 112128 : xlp = atan2(fz,sqrt(x2py2))
835 112128 : slp = sin(xlp)
836 112128 : clp = cos(xlp)
837 112128 : qdlat = xlp*rtd ! output
838 112128 : clm = sqrror*clp ! output
839 112128 : 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 112128 : xlatm = rtd*acos(clm)
845 : !
846 : ! If southern magnetic hemisphere, reverse sign of xlatm
847 : !
848 112128 : if (slp < 0._r8) xlatm = -xlatm
849 448512 : do i=1,3
850 336384 : grclp = cpm*gradx(i) + spm*grady(i)
851 336384 : rgrlp(i) = r*(clp*gradz(i) - slp*grclp)
852 336384 : grclm(i) = sqrror*grclp
853 448512 : clmgrp(i) = sqrror*(cpm*grady(i)-spm*gradx(i))
854 : enddo
855 112128 : grclm(3) = grclm(3) - sqrror*clp/(2._r8*r)
856 :
857 112128 : end subroutine gradlpv
858 : !-----------------------------------------------------------------------
859 112128 : 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 112128 : rr = re + hr
897 112128 : simoslm = 2._r8/sqrt(4._r8 - 3._r8*clm*clm)
898 112128 : sim = simoslm*sin(xlatm*dtr)
899 112128 : bmag = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
900 112128 : d1db = 0._r8
901 112128 : d2db = 0._r8
902 448512 : do i=1,3
903 336384 : bhat(i) = b(i)/bmag
904 336384 : d1(i) = rr*clmgrp(i)
905 336384 : d1db = d1db + d1(i)*bhat(i)
906 336384 : d2(i) = rr*simoslm*grclm(i)
907 448512 : d2db = d2db + d2(i)*bhat(i)
908 : enddo
909 : !
910 : ! Ensure that d1,d2 are exactly perpendicular to B:
911 : !
912 448512 : do i=1,3
913 336384 : d1(i) = d1(i) - d1db*bhat(i)
914 448512 : d2(i) = d2(i) - d2db*bhat(i)
915 : enddo
916 112128 : e3(1) = d1(2)*d2(3) - d1(3)*d2(2)
917 112128 : e3(2) = d1(3)*d2(1) - d1(1)*d2(3)
918 112128 : e3(3) = d1(1)*d2(2) - d1(2)*d2(1)
919 112128 : d = bhat(1)*e3(1) + bhat(2)*e3(2) + bhat(3)*e3(3)
920 448512 : do i=1,3
921 336384 : d3(i) = bhat(i)/d
922 448512 : e3(i) = bhat(i)*d ! Ensure that e3 lies along bhat.
923 : enddo
924 112128 : e1(1) = d2(2)*d3(3) - d2(3)*d3(2)
925 112128 : e1(2) = d2(3)*d3(1) - d2(1)*d3(3)
926 112128 : e1(3) = d2(1)*d3(2) - d2(2)*d3(1)
927 112128 : e2(1) = d3(2)*d1(3) - d3(3)*d1(2)
928 112128 : e2(2) = d3(3)*d1(1) - d3(1)*d1(3)
929 112128 : e2(3) = d3(1)*d1(2) - d3(2)*d1(1)
930 112128 : w = rr*rr*clm*abs(sim)/(bmag*d)
931 112128 : si = -bhat(3)
932 112128 : f1(1) = rgrlp(2)
933 112128 : f1(2) = -rgrlp(1)
934 112128 : f2(1) = -d1(2)*r3_2
935 112128 : f2(2) = d1(1)*r3_2
936 112128 : f = f1(1)*f2(2) - f1(2)*f2(1)
937 :
938 112128 : end subroutine basevec
939 : !-----------------------------------------------------------------------
940 198518784 : 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 198518784 : gpl = sqrt( gb(2 )**2+ gb(3 )**2+ gb(4 )**2)
956 198518784 : ctp = gb(2 )/gpl
957 :
958 198518784 : colat = (acos(ctp))*rtd
959 198518784 : 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 198518784 : vp = .2_r8*gpl*re
965 :
966 198518784 : end subroutine apex_dypol
967 : !-----------------------------------------------------------------------
968 595547136 : 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 198515712 : call cofrm(date)
980 198515712 : call apex_dypol(clatp,polon,vpol)
981 : !
982 : ! colat,ctp,stp,elon,vp are in module data.
983 : !
984 198515712 : colat = clatp
985 198515712 : ctp = cos(clatp*dtr)
986 198515712 : stp = sqrt(1._r8-ctp*ctp)
987 :
988 198515712 : elon = polon
989 198515712 : vp = vpol
990 :
991 198515712 : vmp = 0._r8
992 : !
993 : ! Last 7 args of linapx are output:
994 : !
995 198515712 : call linapx(dlat,dlon,alt, aht,alat,alon,xmag,ymag,zmag,bmag)
996 :
997 198515712 : xmag = xmag*1.e5_r8
998 198515712 : ymag = ymag*1.e5_r8
999 198515712 : zmag = zmag*1.e5_r8
1000 198515712 : bmag = bmag*1.e5_r8
1001 198515712 : call gd2cart (dlat,dlon,alt,x,y,z)
1002 198515712 : iflag = 3
1003 198515712 : xre = x/re ; yre = y/re ; zre = z/re
1004 198515712 : call feldg(iflag,xre,yre,zre,bx,by,bz,vmp)
1005 :
1006 198515712 : end subroutine apex_sub
1007 : !-----------------------------------------------------------------------
1008 11800602624 : 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 198515712 : iflag = 2 ! gclat,r are returned
1039 198515712 : call convrt(iflag,gdlat,alt,gclat,r)
1040 :
1041 198515712 : singml = ctp*sin(gclat*dtr) + stp*cos(gclat*dtr)*cos((glon-elon)*dtr)
1042 198515712 : cgml2 = max(0.25_r8,1._r8-singml*singml)
1043 198515712 : ds = .06_r8*r/cgml2 - 370._r8 ! ds is in module data
1044 :
1045 198515712 : yapx = 0._r8 ! init (module data)
1046 : !
1047 : ! Convert from geodetic to earth centered cartesian coordinates:
1048 : !
1049 198515712 : call gd2cart(gdlat,glon,alt,y(1),y(2),y(3))
1050 198515712 : nstp = 0
1051 : !
1052 : ! Get magnetic field components to determine the direction for tracing field line:
1053 : !
1054 198515712 : iflag = 1
1055 198515712 : call feldg(iflag,gdlat,glon,alt,xmag,ymag,zmag,fmag)
1056 :
1057 198515712 : 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 11616880128 : iflag = 2 ! module data bx,by,bz,bb are returned
1064 11616880128 : y1 = y(1)/re ; y2 = y(2)/re ; y3 = y(3)/re
1065 11616880128 : call feldg(iflag,y1,y2,y3,bx,by,bz,bb)
1066 11616880128 : nstp = nstp + 1
1067 : !
1068 : ! Quit if too many steps.
1069 : !
1070 11616880128 : if (nstp >= maxs) then
1071 14793216 : rho = sqrt(y(1)*y(1) + y(2)*y(2))
1072 14793216 : iflag = 3 ! xlat and ht are returned
1073 14793216 : call convrt(iflag,xlat,ht,rho,y(3))
1074 14793216 : xlon = rtd*atan2(y(2),y(1))
1075 14793216 : iflag = 1
1076 14793216 : call feldg(iflag,xlat,xlon,ht,bnrth,beast,bdown,babs)
1077 14793216 : call dipapx(xlat,xlon,ht,bnrth,beast,bdown,aht,alon)
1078 14793216 : alat = -sgn*rtd*acos(sqrt(1._r8/aht))
1079 14793216 : return
1080 : endif
1081 : !
1082 : ! Find next point using adams algorithm after 7 points
1083 : !
1084 11602086912 : call itrace(iapx)
1085 11602086912 : if (iapx == 1) goto 100
1086 : !
1087 : ! Maximum radius just passed. Find apex coords
1088 : !
1089 183722496 : call fndapx(alt,zmag,aht,alat,alon)
1090 :
1091 : end subroutine linapx
1092 : !-----------------------------------------------------------------------
1093 1925984256 : 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 1925984256 : if (iflag < 3) then ! geodetic to geocentric
1150 : !
1151 : ! Compute rho,z
1152 1360023552 : sinlat = sin(gdlat*dtr)
1153 1360023552 : coslat = sqrt(1._r8-sinlat*sinlat)
1154 1360023552 : d = sqrt(1._r8-e2*sinlat*sinlat)
1155 1360023552 : z = (alt+ome2req/d)*sinlat
1156 1360023552 : rho = (alt+req/d)*coslat
1157 1360023552 : x1 = rho
1158 1360023552 : x2 = z
1159 1360023552 : if (iflag == 1) return
1160 : !
1161 : ! Compute gclat,rkm
1162 198515712 : rkm = sqrt(z*z+rho*rho)
1163 198515712 : gclat = rtd*atan2(z,rho)
1164 198515712 : x1 = gclat
1165 198515712 : x2 = rkm
1166 198515712 : return ! iflag == 2
1167 : endif ! iflag < 3
1168 : !
1169 : ! Geocentric to geodetic:
1170 565960704 : if (iflag == 3) then
1171 565960704 : rho = x1
1172 565960704 : z = x2
1173 565960704 : rkm = sqrt(z*z+rho*rho)
1174 565960704 : scl = z/rkm
1175 565960704 : 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 565960704 : ri = req/rkm
1187 565960704 : a2 = ri*(a21+ri*(a22+ri* a23))
1188 565960704 : a4 = ri*(a41+ri*(a42+ri*(a43+ri*a44)))
1189 565960704 : a6 = ri*(a61+ri*(a62+ri* a63))
1190 565960704 : a8 = ri*(a81+ri*(a82+ri*(a83+ri*a84)))
1191 565960704 : ccl = sqrt(1._r8-scl*scl)
1192 565960704 : s2cl = 2._r8*scl*ccL
1193 565960704 : c2cl = 2._r8*ccl*ccl-1._r8
1194 565960704 : s4cl = 2._r8*s2cl*c2cl
1195 565960704 : c4cl = 2._r8*c2cl*c2cl-1._r8
1196 565960704 : s8cl = 2._r8*s4cl*c4cl
1197 565960704 : s6cl = s2cl*c4cl+c2cl*s4cl
1198 565960704 : dltcl = s2cl*a2+s4cl*a4+s6cl*a6+s8cl*a8
1199 565960704 : gdlat = dltcl*rtd+gclat
1200 565960704 : sgl = sin(gdlat*dtr)
1201 565960704 : alt = rkm*cos(dltcl)-req*sqrt(1._r8-e2*sgl*sgl)
1202 :
1203 : end subroutine convrt
1204 : !-----------------------------------------------------------------------
1205 1161507840 : 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 1161507840 : iflag = 1 ! Convert from geodetic to cylindrical (rho,z are output)
1217 1161507840 : call convrt(iflag,gdlat,alt,rho,z)
1218 :
1219 1161507840 : ang = glon*dtr
1220 1161507840 : x = rho*cos(ang)
1221 1161507840 : y = rho*sin(ang)
1222 :
1223 1161507840 : end subroutine gd2cart
1224 : !-----------------------------------------------------------------------
1225 12579872256 : 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 12579872256 : if (iflag == 1) then
1278 764476416 : is = 1
1279 764476416 : rlat = glat*dtr
1280 764476416 : ct = sin(rlat)
1281 764476416 : st = cos(rlat)
1282 764476416 : rlon = glon*dtr
1283 764476416 : cp = cos(rlon)
1284 764476416 : sp = sin(rlon)
1285 764476416 : call gd2cart(glat,glon,alt,xxx,yyy,zzz)
1286 764476416 : xxx = xxx/re
1287 764476416 : yyy = yyy/re
1288 764476416 : zzz = zzz/re
1289 : else
1290 11815395840 : is = 2
1291 11815395840 : xxx = glat
1292 11815395840 : yyy = glon
1293 11815395840 : zzz = alt
1294 : endif
1295 12579872256 : rq = 1._r8/(xxx**2+yyy**2+zzz**2)
1296 12579872256 : xi(1) = xxx*rq
1297 12579872256 : xi(2) = yyy*rq
1298 12579872256 : xi(3) = zzz*rq
1299 12579872256 : ihmax = nmax*nmax+1
1300 12579872256 : last = ihmax+nmax+nmax
1301 12579872256 : 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 12579872256 : if (iflag /= 3) then
1307 >24391*10^8 : do i=1,last
1308 >24391*10^8 : g(i) = gb(i) ! gb is module data from cofrm
1309 : enddo
1310 : else
1311 39107595264 : do i=1,last
1312 39107595264 : g(i) = gv(i) ! gv is module data from cofrm
1313 : enddo
1314 : endif
1315 :
1316 >35223*10^7 : do i=ihmax,last
1317 >35223*10^7 : h(i) = g(i)
1318 : enddo
1319 :
1320 : mk = 3
1321 : if (imax == 1) mk = 1
1322 :
1323 25159744512 : do k=1,mk,2
1324 25159744512 : i = imax
1325 25159744512 : ih = ihmax
1326 :
1327 : 100 continue
1328 >31449*10^7 : il = ih-i
1329 >31449*10^7 : f = 2._r8/dble(i-k+2)
1330 >31449*10^7 : x = xi(1)*f
1331 >31449*10^7 : y = xi(2)*f
1332 >31449*10^7 : z = xi(3)*(f+f)
1333 :
1334 >31449*10^7 : i = i-2
1335 >31449*10^7 : if (i < 1) then
1336 12579872256 : h(il) = g(il) + z*h(ih) + 2._r8*(x*h(ih+1)+y*h(ih+2))
1337 >30191*10^7 : elseif (i == 1) then
1338 25159744512 : h(il+2) = g(il+2) + z*h(ih+2) + x*h(ih+4) - y*(h(ih+3)+h(ih))
1339 25159744512 : h(il+1) = g(il+1) + z*h(ih+1) + y*h(ih+4) + x*(h(ih+3)-h(ih))
1340 25159744512 : h(il) = g(il) + z*h(ih) + 2._r8*(x*h(ih+1)+y*h(ih+2))
1341 : else
1342 >27675*10^7 : do m=3,i,2
1343 >16605*10^8 : ihm = ih+m
1344 >16605*10^8 : ilm = il+m
1345 >66421*10^8 : h(ilm+1) = g(ilm+1)+ z*h(ihm+1) + x*(h(ihm+3)-h(ihm-1))- &
1346 >83027*10^8 : y*(h(ihm+2)+h(ihm-2))
1347 >33210*10^8 : h(ilm) = g(ilm) + z*h(ihm) + x*(h(ihm+2)-h(ihm-2))+ &
1348 >49816*10^8 : y*(h(ihm+3)+h(ihm-1))
1349 : enddo
1350 >27675*10^7 : h(il+2) = g(il+2) + z*h(ih+2) + x*h(ih+4) - y*(h(ih+3)+h(ih))
1351 >27675*10^7 : h(il+1) = g(il+1) + z*h(ih+1) + y*h(ih+4) + x*(h(ih+3)-h(ih))
1352 >27675*10^7 : h(il) = g(il) + z*h(ih) + 2._r8*(x*h(ih+1)+y*h(ih+2))
1353 : endif
1354 :
1355 >31449*10^7 : ih = il
1356 >32707*10^7 : if (i >= k) goto 100
1357 : enddo ! k=1,mk,2
1358 :
1359 12579872256 : s = .5_r8*h(1)+2._r8*(h(2)*xi(3)+h(3)*xi(1)+h(4)*xi(2))
1360 12579872256 : t = (rq+rq)*sqrt(rq)
1361 12579872256 : bxxx = t*(h(3)-s*xxx)
1362 12579872256 : byyy = t*(h(4)-s*yyy)
1363 12579872256 : bzzz = t*(h(2)-s*zzz)
1364 12579872256 : babs = sqrt(bxxx**2+byyy**2+bzzz**2)
1365 12579872256 : if (is .eq. 1) then ! (convert back to geodetic)
1366 764476416 : beast = byyy*cp-bxxx*sp
1367 764476416 : brho = byyy*sp+bxxx*cp
1368 764476416 : bnrth = bzzz*st-brho*ct
1369 764476416 : bdown = -bzzz*ct-brho*st
1370 11815395840 : elseif (is .eq. 2) then ! (leave in earth centered cartesian)
1371 11815395840 : bnrth = bxxx
1372 11815395840 : beast = byyy
1373 11815395840 : 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 12579872256 : if (iflag == 3) babs = (bxxx*xxx + byyy*yyy + bzzz*zzz)*re*.1_r8
1382 :
1383 12579872256 : end subroutine feldg
1384 : !-----------------------------------------------------------------------
1385 14793216 : 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 14793216 : bhor = sqrt(bnorth*bnorth + beast*beast)
1432 14793216 : if (bhor == 0._r8) then
1433 0 : alon = 0._r8
1434 0 : a = 1.e34_r8
1435 0 : return
1436 : endif
1437 :
1438 14793216 : cottd = bdown*.5_r8/bhor
1439 14793216 : std = 1._r8/sqrt(1._r8+cottd*cottd)
1440 14793216 : ctd = cottd*std
1441 14793216 : sb = -beast/bhor
1442 14793216 : cb = -bnorth/bhor
1443 14793216 : ctg = sin(gdlat*dtr)
1444 14793216 : stg = cos(gdlat*dtr)
1445 14793216 : ang = (gdlon-elon)*dtr
1446 14793216 : sang = sin(ang)
1447 14793216 : cang = cos(ang)
1448 14793216 : cte = ctg*std + stg*ctd*cb
1449 14793216 : ste = sqrt(1._r8 - cte*cte)
1450 14793216 : sa = sb*ctd/ste
1451 14793216 : ca = (std*stg - ctd*ctg*cb)/ste
1452 14793216 : capang = ca*cang - sa*sang
1453 14793216 : sapang = ca*sang + sa*cang
1454 14793216 : stfcpa = ste*ctp*capang - cte*stp
1455 14793216 : stfspa = sapang*ste
1456 14793216 : alon = atan2(stfspa,stfcpa)*rtd
1457 14793216 : r = alt + re
1458 14793216 : ha = alt + r*cottd*cottd
1459 14793216 : a = 1._r8 + ha/req
1460 :
1461 : end subroutine dipapx
1462 : !-----------------------------------------------------------------------
1463 11602086912 : 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 11602086912 : 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 11602086912 : yploc(1,4) = sgn*bx/bb
1487 11602086912 : yploc(2,4) = sgn*by/bb
1488 11602086912 : yploc(3,4) = sgn*bz/bb
1489 :
1490 11602086912 : if (nstp > 7) then
1491 40853084160 : do i=1,3
1492 30639813120 : yapx(i,1) = yapx(i,2)
1493 30639813120 : yapx(i,2) = y(i)
1494 30639813120 : yp(i) = y(i)
1495 30639813120 : term = 55._r8*yploc(i,4)-59._r8*yploc(i,3)+37._r8*yploc(i,2)-9._r8*yploc(i,1)
1496 30639813120 : y(i) = yp(i) + d24*term
1497 30639813120 : yapx(i,3) = y(i)
1498 >13277*10^7 : do j=1,3
1499 >12255*10^7 : yploc(i,j) = yploc(i,j+1)
1500 : enddo
1501 : enddo
1502 10213271040 : rc = rdus ( y(1), y(2), y(3))
1503 10213271040 : rp = rdus (yp(1), yp(2), yp(3))
1504 10213271040 : if (rc < rp) iapx=2
1505 10213271040 : return
1506 : endif
1507 :
1508 5555263488 : do i=1,3
1509 1388815872 : select case (nstp)
1510 : case (1)
1511 595547136 : d2 = ds/2._r8
1512 595547136 : d6 = ds/6._r8
1513 595547136 : d12 = ds/12._r8
1514 595547136 : d24 = ds/24._r8
1515 595547136 : yploc(i,1)= yploc(i,4)
1516 595547136 : yp(i) = y(i)
1517 595547136 : yapx(i,1) = y(i)
1518 595547136 : y(i) = yp(i) + ds*yploc(i,1)
1519 : case (2)
1520 595547136 : yploc(i,2) = yploc(i,4)
1521 595547136 : y(i) = yp(i) + d2*(yploc(i,2)+yploc(i,1))
1522 : case (3)
1523 595547136 : y(i) = yp(i) + d6*(2._r8*yploc(i,4)+yploc(i,2)+3._r8*yploc(i,1))
1524 : case (4)
1525 595547136 : yploc(i,2) = yploc(i,4)
1526 595547136 : yapx(i,2) = y(i)
1527 595547136 : yp(i) = y(i)
1528 595547136 : y(i) = yp(i) + d2*(3._r8*yploc(i,2)-yploc(i,1))
1529 : case (5)
1530 595547136 : y(i) = yp(i) + d12*(5._r8*yploc(i,4)+8._r8*yploc(i,2)-yploc(i,1))
1531 : case (6)
1532 595547136 : yploc(i,3) = yploc(i,4)
1533 595547136 : yp(i) = y(i)
1534 595547136 : yapx(i,3) = y(i)
1535 595547136 : y(i) = yp(i) + d12*(23._r8*yploc(i,3)-16._r8*yploc(i,2)+5._r8*yploc(i,1))
1536 : case (7)
1537 593164800 : yapx(i,1) = yapx(i,2)
1538 593164800 : yapx(i,2) = yapx(i,3)
1539 593164800 : y(i) = yp(i) + d24*(9._r8*yploc(i,4)+19._r8*yploc(i,3)-5._r8*yploc(i,2)+yploc(i,1))
1540 593164800 : yapx(i,3) = y(i)
1541 : case default
1542 0 : write(iulog,"('>>> itrace: unresolved case nstp=',i4)") nstp
1543 4166447616 : call endrun( 'itrace' )
1544 : end select
1545 : enddo
1546 : !
1547 : ! Signal if apex passed:
1548 : !
1549 1388815872 : if (nstp == 6 .or. nstp == 7) then
1550 396237312 : rc = rdus( yapx(1,3), yapx(2,3), yapx(3,3))
1551 396237312 : rp = rdus( yapx(1,2), yapx(2,2), yapx(3,2))
1552 396237312 : if (rc < rp) iapx=2
1553 : endif
1554 :
1555 : end subroutine itrace
1556 : !-----------------------------------------------------------------------
1557 21219016704 : real(r8) function rdus(d,e,f)
1558 : real(r8),intent(in) :: d,e,f
1559 21219016704 : rdus = sqrt(d**2 + e**2 + f**2)
1560 21219016704 : end function rdus
1561 : !-----------------------------------------------------------------------
1562 183722496 : 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 183722496 : iflag_feldg = 1
1578 183722496 : iflag_convrt = 3
1579 734889984 : do i=1,3
1580 551167488 : rho = sqrt(yapx(1,i)**2+yapx(2,i)**2)
1581 551167488 : call convrt(iflag_convrt,gdlt,ht(i),rho,yapx(3,i))
1582 551167488 : gdln = rtd*atan2(yapx(2,i),yapx(1,i))
1583 734889984 : 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 734889984 : do i=1,3
1589 734889984 : 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 183722496 : 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 183722496 : xinter = max(alt,xinter)
1598 183722496 : 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 183722496 : if (a < 1._r8) then
1604 0 : write(iulog,"('>>> fndapx: a=',e12.4,' < 1.')") a
1605 0 : call endrun( 'fndapx' )
1606 : endif
1607 :
1608 183722496 : rasq = rtd*acos(sqrt(1._r8/a))
1609 183722496 : 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 183722496 : xlon = atan2(yloc(2),yloc(1))
1627 183722496 : ang = xlon-elon*dtr
1628 183722496 : cang = cos(ang)
1629 183722496 : sang = sin(ang)
1630 183722496 : r = sqrt(yloc(1)**2+yloc(2)**2+yloc(3)**2)
1631 183722496 : cte = yloc(3)/r
1632 183722496 : ste = sqrt(1._r8-cte*cte)
1633 183722496 : stfcpa = ste*ctp*cang - cte*stp
1634 183722496 : stfspa = sang*ste
1635 183722496 : alon = atan2(stfspa,stfcpa)*rtd
1636 :
1637 183722496 : end subroutine fndapx
1638 : !-----------------------------------------------------------------------
1639 734889984 : 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 734889984 : (a1-a2)*(a7-a1)*(a7-a2)*a6)/((a1-a2)*(a1-a3)*(a2-a3))
1649 :
1650 734889984 : 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 113280 : 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 113280 : ier = 0
1705 113280 : glonloc = glon
1706 113280 : if (glonloc < gplon(1)) glonloc = glonloc + 360._r8
1707 113280 : if (glonloc > gplon(nlon)) glonloc = glonloc - 360._r8
1708 : !
1709 10251072 : i0 = 0
1710 10251072 : do i=1,nlat-1
1711 10251072 : if (glat >= gplat(i).and.glat <= gplat(i+1)) then
1712 113280 : i0 = i
1713 113280 : dlat = gplat(i+1)-gplat(i)
1714 113280 : xi = (glat - gplat(i)) / dlat
1715 113280 : exit
1716 : endif
1717 : enddo
1718 113280 : 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 20505248 : j0 = 0
1726 20505248 : do j=1,nlon-1
1727 20505248 : if (glon >= gplon(j).and.glon <= gplon(j+1)) then
1728 113280 : j0 = j
1729 113280 : dlon = gplon(j+1)-gplon(j)
1730 113280 : yj = (glon - gplon(j)) / dlon
1731 113280 : exit
1732 : endif
1733 : enddo
1734 113280 : 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 113280 : k0 = 0
1742 113280 : do k=1,nalt-1
1743 113280 : if (alt >= gpalt(k).and.alt <= gpalt(k+1)) then
1744 113280 : k0 = k
1745 113280 : hti = re/(re+alt)
1746 113280 : diht = re/(re+gpalt(k+1)) - re/(re+gpalt(k))
1747 113280 : zk = (hti - re/(re+gpalt(k))) / diht
1748 113280 : exit
1749 : endif
1750 : enddo
1751 113280 : 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 1699200 : call trilin(xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fx,dfxdn,dfxde,dfxdd)
1759 113280 : dfxdth = -dfxdn*rtd/dlat
1760 113280 : dfxdln = dfxde*rtd/dlon
1761 113280 : dfxdh = -hti*hti*dfxdd/(re*diht)
1762 :
1763 1699200 : call trilin(yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fy,dfydn,dfyde,dfydd)
1764 113280 : dfydth = -dfydn*rtd/dlat
1765 113280 : dfydln = dfyde*rtd/dlon
1766 113280 : dfydh = -hti*hti*dfydd/(re*diht)
1767 :
1768 1699200 : call trilin(zarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fz,dfzdn,dfzde,dfzdd)
1769 113280 : dfzdth = -dfzdn*rtd/dlat
1770 113280 : dfzdln = dfzde*rtd/dlon
1771 113280 : dfzdh = -hti*hti*dfzdd/(re*diht)
1772 :
1773 1699200 : call trilin(varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,fv,dfvdn,dfvde,dfvdd)
1774 113280 : dfvdth = -dfvdn*rtd/dlat
1775 113280 : dfvdln = dfvde*rtd/dlon
1776 113280 : dfvdh = -hti*hti*dfvdd/(re*diht)
1777 :
1778 113280 : if (nlat < 3) return
1779 : !
1780 : ! Improve calculation of longitudinal derivatives near poles
1781 : !
1782 113280 : if (glat < dlat-90._r8) then
1783 1728 : fac = .5_r8*xi
1784 1728 : omfac = 1._r8 - fac
1785 1728 : xi = xi - 1._r8
1786 1728 : i0 = i0 + 1
1787 25920 : call trilin (xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1788 1728 : dfxdln = dfxdln*omfac + fac*dmdfde*rtd/dlon
1789 25920 : call trilin (yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1790 1728 : dfydln = dfydln*omfac + fac*dmdfde*rtd/dlon
1791 25920 : call trilin (varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1792 1728 : dfvdln = dfvdln*omfac + fac*dmdfde*rtd/dlon
1793 : endif
1794 :
1795 113280 : if (glat > 90._r8-dlat) then
1796 1728 : fac = .5_r8*(1._r8-xi)
1797 1728 : omfac = 1._r8 - fac
1798 1728 : xi = xi + 1._r8
1799 1728 : i0 = i0 - 1
1800 25920 : call trilin (xarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1801 1728 : dfxdln = dfxdln*omfac + fac*dmdfde*rtd/dlon
1802 25920 : call trilin (yarray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1803 1728 : dfydln = dfydln*omfac + fac*dmdfde*rtd/dlon
1804 25920 : call trilin (varray(i0:i0+1,j0:j0+1,k0:k0+1),xi,yj,zk,dmf,dmdfdn,dmdfde,dmdfdd)
1805 1728 : dfvdln = dfvdln*omfac + fac*dmdfde*rtd/dlon
1806 : endif
1807 :
1808 : end subroutine intrp
1809 : !-----------------------------------------------------------------------
1810 463488 : 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 463488 : omxi = 1._r8 - xi
1834 463488 : omyj = 1._r8 - yj
1835 463488 : 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 463488 : + 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 463488 : + (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 463488 : + (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 463488 : + (u(2,2,2)-u(2,2,1))*xi*yj
1858 :
1859 463488 : end subroutine trilin
1860 : !-----------------------------------------------------------------------
1861 113280 : 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 113280 : cph = cos((glon-elon)*dtr)
1878 113280 : sph = sin((glon-elon)*dtr)
1879 113280 : cth = sin(glat*dtr)
1880 113280 : sth = cos(glat*dtr)
1881 113280 : ctm = ctp*cth + stp*sth*cph
1882 113280 : fx = fx + sth*ctp*cph - cth*stp
1883 113280 : fy = fy + sth*sph
1884 113280 : fz = fz + ctm
1885 113280 : fv = fv - ctm
1886 :
1887 113280 : dfxdth = dfxdth + ctp*cth*cph + stp*sth
1888 113280 : dfydth = dfydth + cth*sph
1889 113280 : dfzdth = dfzdth - ctp*sth + stp*cth*cph
1890 113280 : dfvdth = dfvdth + ctp*sth - stp*cth*cph
1891 :
1892 113280 : dfxdln = dfxdln - ctp*sth*sph
1893 113280 : dfydln = dfydln + sth*cph
1894 113280 : dfzdln = dfzdln - stp*sth*sph
1895 113280 : dfvdln = dfvdln + stp*sth*sph
1896 :
1897 113280 : 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 198517248 : 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 397034496 : real(r8) :: gh(n1*ncn1 + n2*ncn2 + 1)
1943 :
1944 : real(r8),parameter :: alt = 0._r8
1945 : integer, parameter :: isv=0
1946 :
1947 198517248 : ngh = n1*ncn1 + n2*ncn2 + 1 ! not sure why the extra +1
1948 :
1949 198517248 : 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 198517248 : 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 3970344960 : do n=1,ncn1
1966 3771827712 : i = (n-1)*n1
1967 >45658*10^7 : 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 1389620736 : do n=1,ncn2
1971 1191103488 : i = n1*ncn1 + (n-1)*n2
1972 >23365*10^7 : 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 198517248 : gh(ngh) = 0._r8 ! not sure why gh is dimensioned with the extra element, so set it to 0.
1976 :
1977 198517248 : if (date < apex_end_yr-10) then
1978 198517248 : t = 0.2_r8*(date - year1)
1979 198517248 : ll = t
1980 198517248 : one = ll
1981 198517248 : t = t - one
1982 198517248 : if (date < year2-5) then
1983 198517248 : nmx = 10
1984 198517248 : nc = nmx*(nmx+2)
1985 198517248 : ll = nc*ll
1986 198517248 : 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 198517248 : 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 198517248 : l = 1
2012 198517248 : m = 1
2013 198517248 : n = 0
2014 : !
2015 : ! Set outputs gb(ncoef) and gv(ncoef)
2016 : ! These are module data above.
2017 : !
2018 198517248 : gb(:) = 0._r8
2019 198517248 : gv(:) = 0._r8
2020 198517248 : f0 = -1.e-5_r8
2021 13102138368 : do k=2,kmx
2022 12903621120 : if (n < m) then
2023 1985172480 : m = 0
2024 1985172480 : n = n+1
2025 : endif ! n < m
2026 12903621120 : lm = ll + l
2027 12903621120 : if (m == 0) f0 = f0 * dble(n)/2._r8
2028 1985172480 : if (m == 0) f = f0 / sqrt(2.0_r8)
2029 12903621120 : nn = n+1
2030 :
2031 12903621120 : if (m /= 0) then
2032 10918448640 : f = f / sqrt(dble(n-m+1) / dble(n+m) )
2033 10918448640 : gb(l+1) = (tc*gh(lm) + t*gh(lm+nc))* f
2034 : else
2035 1985172480 : gb(l+1) = (tc*gh(lm) + t*gh(lm+nc))* f0
2036 : endif
2037 12903621120 : gv(l+1) = gb(l+1)/dble(nn)
2038 12903621120 : if (m /= 0) then
2039 10918448640 : gb(l+2) = (tc*gh(lm+1) + t*gh(lm+nc+1))*f
2040 10918448640 : gv(l+2) = gb(l+2)/dble(nn)
2041 10918448640 : l = l+2
2042 : else
2043 : l = l+1
2044 : endif
2045 13102138368 : 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 198517248 : end subroutine cofrm
2052 : !-----------------------------------------------------------------------
2053 72960 : 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 72960 : sbsllat=0._r8 ; sbsllon=0._r8
2095 :
2096 72960 : 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 72960 : nleap = (iyr-1601)/4
2101 72960 : nleap = nleap - 99
2102 72960 : 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 72960 : 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 72960 : 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 72960 : g0 = -2.472_r8 + (-.2558905_r8*(yr-4*nleap) - 3.79617e-2_r8*nleap)
2137 : !
2138 : ! Universal time in seconds:
2139 72960 : ut = dble(ihr*3600 + imn*60) + sec
2140 : !
2141 : ! Days (including fraction) since 12 UT on January 1 of IYR:
2142 72960 : df = (ut/86400._r8 - 1.5_r8) + iday
2143 : !
2144 : ! Addition to Mean longitude of Sun since January 1 of IYR:
2145 72960 : lf = .9856474_r8*df
2146 : !
2147 : ! Addition to Mean anomaly since January 1 of IYR:
2148 72960 : gf = .9856003_r8*df
2149 : !
2150 : ! Mean longitude of Sun:
2151 72960 : l = l0 + lf
2152 : !
2153 : ! Mean anomaly:
2154 72960 : g = g0 + gf
2155 72960 : grad = g*d2r
2156 : !
2157 : ! Ecliptic longitude:
2158 72960 : lambda = l + 1.915_r8*sin(grad) + .020_r8*sin(2._r8*grad)
2159 72960 : lamrad = lambda*d2r
2160 72960 : sinlam = sin(lamrad)
2161 : !
2162 : ! Days (including fraction) since 12 UT on January 1 of 2000:
2163 72960 : n = df + 365._r8*yr + dble(nleap)
2164 : !
2165 : ! Obliquity of ecliptic:
2166 72960 : epsilon = 23.439_r8 - 4.e-7_r8*n
2167 72960 : epsrad = epsilon*d2r
2168 : !
2169 : ! Right ascension:
2170 72960 : alpha = atan2(cos(epsrad)*sinlam,cos(lamrad))*r2d
2171 : !
2172 : ! Declination:
2173 72960 : delta = asin(sin(epsrad)*sinlam)*r2d
2174 : !
2175 : ! Subsolar latitude (output argument):
2176 72960 : sbsllat = delta
2177 : !
2178 : ! Equation of time (degrees):
2179 72960 : etdeg = l - alpha
2180 72960 : nrot = nint(etdeg/360._r8)
2181 72960 : etdeg = etdeg - dble(360*nrot)
2182 : !
2183 : ! Apparent time (degrees):
2184 : ! Earth rotates one degree every 240 s.
2185 72960 : aptime = ut/240._r8 + etdeg
2186 : !
2187 : ! Subsolar longitude (output argument):
2188 72960 : sbsllon = 180._r8 - aptime
2189 72960 : nrot = nint(sbsllon/360._r8)
2190 72960 : sbsllon = sbsllon - dble(360*nrot)
2191 :
2192 72960 : end subroutine apex_subsol
2193 :
2194 : !-----------------------------------------------------------------------
2195 1050624 : 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 1050624 : ctp = cos(colat*dtor)
2229 1050624 : stp = sqrt(1._r8 - ctp*ctp)
2230 1050624 : ang = (xlon-elon)*dtor
2231 1050624 : cang = cos(ang)
2232 1050624 : sang = sin(ang)
2233 1050624 : cte = sin(xlat*dtor)
2234 1050624 : ste = sqrt(1._r8-cte*cte)
2235 1050624 : stfcpa = ste*ctp*cang - cte*stp
2236 1050624 : stfspa = sang*ste
2237 1050624 : mlon = atan2(stfspa,stfcpa)*rtod
2238 :
2239 1050624 : end subroutine solgmlon
2240 : !-----------------------------------------------------------------------
2241 : !-----------------------------------------------------------------------
2242 1050624 : 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 1050624 : call solgmlon (sbsllat,sbsllon,clatp,polon,smlon)
2274 1050624 : mlt = (alon - smlon)/15.0_r8 + 12.0_r8
2275 1050624 : if (mlt .ge. 24.0_r8) mlt = mlt - 24.0_r8
2276 1050624 : if (mlt .lt. 0._r8) mlt = mlt + 24.0_r8
2277 1050624 : return
2278 : end subroutine apex_magloctm
2279 :
2280 : !-----------------------------------------------------------------------
2281 : end module apex
|
