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