Line data Source code
1 : module wei05sc
2 : !
3 : ! The Weimer model of high-latitude potential created by Daniel Weimer and
4 : ! if extracted, distributed, or used for any purpose other than as implemented
5 : ! in the NCAR TIEGCM and CESM/WACCM models, please contact Dan Weimer for
6 : ! further information and discussion.
7 : !
8 : ! 2005 Version of the electric and magnetic potential (FAC) models
9 : ! by Dan Weimer. Uses Spherical Cap Harmonic Analysis (SCHA) functions.
10 : ! Model description is in:
11 : ! Weimer, D. R., Predicting Surface Geomagnetic Variations Using Ionospheric
12 : ! Electrodynamic Models, J. Geophys. Res., 110, A12307, doi:10.1029/
13 : ! 2005JA011270, 2005.
14 : ! Some information about the model (such as outer boundary calculation)
15 : ! is also in the earlier paper:
16 : ! Weimer, D. R. (2005), Improved ionospheric electrodynamic models and
17 : ! application to calculating Joule heating rates, J. Geophys. Res., 110,
18 : ! A05306, doi:10.1029/2004JA010884.
19 : !
20 : ! For information about the SCHA, see the paper:
21 : ! Haines, G. V., Spherical cap harmonic analysis, J. Geophys. Res., 90, B3,
22 : ! 2583, 1985. (Note that this is in JGR-B, "Solid Earth", rather than JGR-A)
23 : !
24 : ! April, 2008:
25 : ! This f90 module of the Electric Potential model was translated
26 : ! from the original IDL by Ben Foster (NCAR, foster@ucar.edu)
27 : ! Netcdf data file wei05sc.nc was written from original IDL save files
28 : ! W05scBndy.xdr, W05scEpot.xdr, W05scBpot.xdr, and SCHAtable.dat
29 : !
30 : ! September, 2015 btf:
31 : ! Modified for free-format fortran, and for CESM/WACCM (r8, etc).
32 : !
33 : use shr_kind_mod, only: r8 => shr_kind_r8
34 : use shr_kind_mod, only: shr_kind_cl
35 : use spmd_utils, only: masterproc
36 : use cam_logfile, only: iulog
37 : use cam_abortutils, only: endrun
38 : use time_manager, only: get_curr_date
39 : use edyn_maggrid, only: nmlat,nmlon,nmlonp1
40 :
41 : use edyn_maggrid, only: &
42 : ylonm, & ! magnetic latitudes (nmlat) (radians)
43 : ylatm ! magnetic longtitudes (nmlonp1) (radians)
44 : use edyn_solve, only: &
45 : nmlat0, & ! (nmlat+1)/2
46 : phihm ! output: high-latitude potential (nmlonp1,nmlat)
47 :
48 : use physconst, only: pi
49 : use aurora_params, only: aurora_params_set, hpower, ctpoten, theta0
50 : use aurora_params, only: offa, dskofa, dskofc, phid, rrad, offc, phin
51 : implicit none
52 : private
53 :
54 : !
55 : ! Coefficients read from netcdf data file wei05sc.nc:
56 : !
57 : integer,parameter :: &
58 : na=6, nb=7, nex=2, n1_scha=19, n2_scha=7, n3_scha=68, &
59 : csize=28, n_schfits=15, n_alschfits=18
60 : integer :: maxk_scha, maxm_scha, maxl_pot, maxm_pot
61 : real(r8) :: bndya(na), bndyb(nb), ex_bndy(nex), ex_epot(nex),ex_bpot(nex)
62 : real(r8) :: th0s(n3_scha), allnkm(n1_scha,n2_scha,n3_scha)
63 : integer :: ab(csize), ls(csize), ms(csize)
64 : real(r8) :: epot_alschfits(n_alschfits,csize)
65 : real(r8) :: bpot_alschfits(n_alschfits,csize)
66 : real(r8) :: epot_schfits(n_schfits,csize)
67 : real(r8) :: bpot_schfits(n_schfits,csize)
68 : !
69 : ! Intermediate calculations:
70 : !
71 : integer,parameter :: mxtablesize=500
72 : real(r8) :: rad2deg,deg2rad ! set by setmodel
73 : real(r8) :: bndyfitr ! calculated by setboundary
74 : real(r8) :: esphc(csize),bsphc(csize) ! calculated by setmodel
75 : real(r8) :: tmat(3,3) ! from setboundary
76 : real(r8) :: plmtable(mxtablesize,csize),colattable(mxtablesize)
77 : real(r8) :: nlms(csize)
78 : real(r8),allocatable :: wei05sc_fac(:,:) ! field-aligned current output
79 :
80 : ! 05/08 bae: Have ctpoten from both hemispheres from Weimer
81 : real(r8) :: weictpoten(2),phimin,phimax
82 :
83 : !
84 : ! Several items in the public list are for efield.F90 (chemistry/mozart)
85 : ! (dpie_coupling calls the weimer05 driver, but efield calls the individual
86 : ! routines, not the driver)
87 : !
88 : public :: weimer05
89 : public :: weimer05_init
90 :
91 : real(r8), parameter :: r2d = 180._r8/pi ! radians to degrees
92 : real(r8), parameter :: d2r = pi/180._r8 ! degrees to radians
93 :
94 : logical :: debug = .false.
95 :
96 : contains
97 :
98 : !-----------------------------------------------------------------------
99 0 : subroutine weimer05_init(wei05_ncfile)
100 : use infnan, only: nan, assignment(=)
101 :
102 : character(len=*),intent(in) :: wei05_ncfile
103 :
104 0 : allocate(wei05sc_fac(nmlonp1,nmlat))
105 :
106 0 : hpower = nan
107 0 : ctpoten = nan
108 0 : phin = nan
109 0 : phid = nan
110 0 : theta0 = nan
111 0 : offa = nan
112 0 : dskofa = nan
113 0 : rrad = nan
114 0 : offc = nan
115 0 : dskofc = nan
116 :
117 0 : bndya = nan
118 0 : bndyb = nan
119 0 : ex_bndy = nan
120 0 : ex_bpot = nan
121 0 : th0s = nan
122 0 : allnkm = nan
123 0 : bpot_schfits = nan
124 0 : bpot_alschfits = nan
125 :
126 0 : if (wei05_ncfile.ne.'NONE') then
127 0 : call read_wei05_ncfile(wei05_ncfile)
128 0 : aurora_params_set = .true.
129 : endif
130 :
131 0 : end subroutine weimer05_init
132 :
133 : !-----------------------------------------------------------------------
134 0 : subroutine weimer05(by, bz_in, swvel, swden, sunlon)
135 : !
136 : ! 9/16/15 btf: Driver to call Weimer 2005 model for waccm[x].
137 : !
138 :
139 : implicit none
140 : !
141 : ! Args:
142 : real(r8), intent(in) :: bz_in, by, swvel, swden
143 : real(r8), intent(in) :: sunlon
144 :
145 : !
146 : ! Local:
147 :
148 : real(r8) :: angl, angle, bt
149 : integer :: i, j
150 : real(r8) :: rmlt, mlat, tilt, htilt, hem, ut, secs
151 : real(r8), parameter :: fill = 0._r8
152 : integer :: iyear, imon, iday, isecs
153 : real(r8) :: bz
154 :
155 0 : bz = bz_in
156 :
157 0 : hpower = hp_from_bz_swvel(bz,swvel)
158 : !
159 : ! Get current date and time:
160 : !
161 0 : call get_curr_date(iyear,imon,iday,isecs)
162 : !
163 : ! Get sun's location (longitude at all latitudes):
164 : !
165 0 : secs = real(isecs, r8)
166 :
167 : !
168 : ! At least one of by,bz must be non-zero:
169 0 : if (by==0._r8 .and. bz==0._r8) then
170 0 : if (masterproc) then
171 : write(iulog,"(/,'>>> WARNING: by and bz cannot both be zero',&
172 0 : ' when calling the Weimer model: am setting bz=0.01')")
173 : end if
174 0 : bz = 0.01_r8
175 : end if
176 : !
177 0 : bt = sqrt(by**2 + bz**2)
178 0 : angl = atan2(by,bz) * r2d
179 : !
180 : ! Convert from day-of-year to month,day and get tilt from date and ut:
181 : !
182 0 : ut = secs / 3600._r8 ! decimal hours
183 : !
184 : ! Given year and day-of-year, cvt2md returns month and day of month.
185 : ! We do not need this, since get_curr_date returns month and day of month.
186 : ! call cvt2md(iulog,iyear,idoy,imon,iday) ! given iyear,idoy, return imo,ida
187 : !
188 0 : if (debug .and. masterproc) then
189 : write(iulog,"('weimer05: iyear,imon,iday=',3i5,' ut=',f8.2)") &
190 0 : iyear,imon,iday,ut
191 : end if
192 0 : tilt = get_tilt(iyear,imon,iday,ut)
193 0 : if (debug .and. masterproc) then
194 0 : write(iulog,"('weimer05: tilt=',e12.4)") tilt
195 : end if
196 :
197 0 : phihm = 0._r8 ! whole-array init (nmlonp1,nmlat)
198 : !
199 : ! Call Weimer model for southern hemisphere electric potential:
200 : !
201 0 : hem = -1._r8
202 0 : htilt = hem * tilt
203 0 : angle = hem * angl
204 0 : if (debug .and. masterproc) then
205 0 : write(iulog,"('weimer05 call setmodel for SH potential')")
206 : end if
207 0 : call setmodel(angle, bt, htilt, swvel, swden, 'epot')
208 0 : if (debug .and. masterproc) then
209 0 : write(iulog,"('weimer05 after setmodel for SH potential')")
210 : end if
211 0 : do j = 1, nmlat0 ! Spole to equator
212 0 : do i = 1, nmlon
213 : !
214 : ! sunlon: sun's longitude in dipole coordinates
215 : !
216 0 : rmlt = (ylonm(i)-sunlon) * r2d / 15._r8 + 12._r8
217 0 : mlat = abs(ylatm(j))*r2d
218 : !
219 : ! Obtain electric potential and convert from kV to V
220 : !
221 0 : call epotval(mlat,rmlt,fill,phihm(i,j))
222 0 : phihm(i,j) = phihm(i,j)*1000._r8
223 : end do ! i=1,nmlon
224 : end do ! j=1,nmlat0
225 0 : if (debug) write(iulog,"('weimer05: SH phihm min,max=',2es12.4)") &
226 0 : minval(phihm(1:nmlon,1:nmlat0)),maxval(phihm(1:nmlon,1:nmlat0))
227 : !
228 : ! Re-calculate SH values of offa, dskofa, arad, and phid and phin from
229 : ! Weimer 2005 setboundary values of offc, dskofc, and theta0
230 : !
231 0 : call wei05loc (1, by, hpower, sunlon)
232 : !
233 : ! Call Weimer model for southern hemisphere fac:
234 : !
235 0 : if (debug .and. masterproc) then
236 0 : write(iulog,"('weimer05 call setmodel for SH fac')")
237 : end if
238 0 : call setmodel(angle,bt,htilt,swvel,swden,'bpot')
239 0 : if (debug .and. masterproc) then
240 0 : write(iulog,"('weimer05 after setmodel for SH fac')")
241 : end if
242 0 : do j = 1, nmlat0
243 0 : do i = 1, nmlon
244 0 : rmlt = (ylonm(i)-sunlon) * r2d / 15._r8 + 12._r8
245 0 : mlat = abs(ylatm(j))*r2d
246 0 : call mpfac(mlat,rmlt,fill,wei05sc_fac(i,j))
247 : end do ! i=1,nmlon
248 : end do ! j=1,nmlat0
249 : !
250 : ! Call Weimer model for northern hemisphere epot:
251 : !
252 0 : hem = 1._r8
253 0 : htilt = hem * tilt
254 0 : angle = hem * angl
255 0 : if (debug .and. masterproc) then
256 0 : write(iulog,"('weimer05 call setmodel for NH potential')")
257 : end if
258 0 : call setmodel(angle,bt,htilt,swvel,swden,'epot')
259 0 : if (debug .and. masterproc) then
260 0 : write(iulog,"('weimer05 after setmodel for NH potential')")
261 : end if
262 0 : do j = nmlat0+1, nmlat
263 0 : do i = 1, nmlon
264 : !
265 : ! sunlon: sun's longitude in dipole coordinates
266 0 : rmlt = ((ylonm(i) - sunlon) * r2d / 15._r8) + 12._r8
267 0 : mlat = abs(ylatm(j)) * r2d
268 : !
269 : ! Obtain electric potential and convert from kV to V
270 0 : call epotval(mlat, rmlt, fill, phihm(i,j))
271 0 : phihm(i,j) = phihm(i,j) * 1000._r8
272 : end do ! i=1,nmlon
273 : end do ! j=1,nmlat0+1,nmlat
274 0 : if (debug .and. masterproc) then
275 : write(iulog,"('weimer05: NH phihm min,max=',2es12.4)") &
276 0 : minval(phihm(1:nmlon,nmlat0+1:nmlat)), &
277 0 : maxval(phihm(1:nmlon,nmlat0+1:nmlat))
278 : end if
279 : !
280 : ! Re-calculate NH values of offa, dskofa, arad, and Heelis phid and phin
281 : ! from Weimer 2005 setboundary values of offc, dskofc, and theta0
282 : !
283 0 : call wei05loc (2, by, hpower, sunlon)
284 : !
285 : ! Call Weimer model for northern hemisphere fac:
286 0 : if (debug .and. masterproc) then
287 0 : write(iulog,"('weimer05 call setmodel for NH fac')")
288 : end if
289 0 : call setmodel(angle,bt,htilt,swvel,swden,'bpot')
290 0 : if (debug .and. masterproc) then
291 0 : write(iulog,"('weimer05 after setmodel for NH fac')")
292 : end if
293 0 : do j = nmlat0+1, nmlat
294 0 : do i = 1, nmlon
295 0 : rmlt = ((ylonm(i)-sunlon) * r2d / 15._r8) + 12._r8
296 0 : mlat = abs(ylatm(j))*r2d
297 0 : call mpfac(mlat,rmlt,fill,wei05sc_fac(i,j))
298 : end do ! i=1,nmlon
299 : end do ! j=1,nmlat0
300 : !
301 : ! Periodic points:
302 0 : do j = 1, nmlat
303 0 : phihm(nmlonp1,j) = phihm(1,j)
304 0 : wei05sc_fac(nmlonp1,j) = wei05sc_fac(1,j)
305 : end do ! j=1,nmlat
306 : !
307 : ! Calculate ctpoten for each hemisphere:
308 : ! South:
309 : !
310 0 : phimax = -1.e36_r8
311 0 : phimin = 1.e36_r8
312 0 : do j = 1, nmlat0 ! SH
313 0 : do i = 1, nmlon
314 0 : if (phihm(i,j) > phimax) phimax = phihm(i,j)
315 0 : if (phihm(i,j) < phimin) phimin = phihm(i,j)
316 : end do
317 : end do
318 0 : weictpoten(1) = 0.001_r8 * (phimax - phimin)
319 : !
320 : ! North:
321 : !
322 0 : phimax = -1.e36_r8
323 0 : phimin = 1.e36_r8
324 0 : do j = nmlat0+1, nmlat ! NH
325 0 : do i = 1, nmlon
326 0 : if (phihm(i,j) > phimax) phimax = phihm(i,j)
327 0 : if (phihm(i,j) < phimin) phimin = phihm(i,j)
328 : end do
329 : end do
330 0 : weictpoten(2) = 0.001_r8 * (phimax - phimin)
331 : !
332 : ! average of the SH and NH in ctpoten
333 0 : ctpoten = 0.5_r8*(weictpoten(1)+weictpoten(2))
334 :
335 0 : if (masterproc) then
336 : write(iulog,"(a,f8.2,a,2es12.4)") &
337 0 : 'weimer05: ctpoten=', ctpoten, ', phihm min,max=', &
338 0 : minval(phihm), maxval(phihm)
339 : end if
340 : !
341 :
342 0 : end subroutine weimer05
343 : !-----------------------------------------------------------------------
344 0 : subroutine read_wei05_ncfile(file)
345 :
346 : use ioFileMod, only: getfil
347 : use cam_pio_utils, only: cam_pio_openfile, cam_pio_closefile
348 : use pio, only: file_desc_t, pio_nowrite, pio_inq_dimid
349 : use pio, only: pio_inquire_dimension, pio_inq_varid, pio_get_var
350 : !
351 : ! Read coefficients and other data from netcdf data file.
352 : !
353 : ! Arg:
354 : character(len=*), intent(in) :: file
355 : !
356 : ! Local:
357 : integer :: istat
358 : integer :: rd_na, rd_nb, rd_nex, rd_n1_scha, rd_n2_scha, rd_n3_scha
359 : integer :: rd_csize, rd_n_schfits, rd_n_alschfits
360 : integer :: id
361 : character(len=shr_kind_cl) :: filen
362 : character(len=shr_kind_cl) :: errmsg
363 : character(len=*), parameter :: prefix = 'read_wei05_ncfile: '
364 : type(file_desc_t) :: ncid
365 : !
366 : ! Open netcdf file for reading:
367 : !
368 0 : call getfil( file, filen, 0 )
369 0 : call cam_pio_openfile(ncid, filen, PIO_NOWRITE)
370 :
371 0 : if (masterproc) then
372 0 : write(iulog,"('wei05sc: opened netcdf data file',a)") trim(filen)
373 : end if
374 : !
375 : ! Read and check dimensions:
376 : !
377 : ! na=6
378 0 : istat = pio_inq_dimid(ncid, 'na', id)
379 0 : istat = pio_inquire_dimension(ncid, id, len=rd_na)
380 0 : if (rd_na /= na) then
381 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_na /= na: rd_na = ', rd_na,' na = ', na
382 0 : write(iulog,*) trim(errmsg)
383 0 : call endrun(errmsg)
384 : end if
385 : !
386 : ! nb=7
387 : !
388 0 : istat = pio_inq_dimid(ncid, 'nb', id)
389 0 : istat = pio_inquire_dimension(ncid, id, len=rd_nb)
390 0 : if (rd_nb /= nb) then
391 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_nb /= nb: rd_nb = ', rd_nb,' nb = ', nb
392 0 : write(iulog,*) trim(errmsg)
393 0 : call endrun(errmsg)
394 : end if
395 : !
396 : ! nex=2
397 : !
398 0 : istat = pio_inq_dimid(ncid, 'nex', id)
399 0 : istat = pio_inquire_dimension(ncid, id, len=rd_nex)
400 0 : if (rd_nex /= nex) then
401 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_nex /= nex rd_nex = ', rd_nex,' nex = ', nex
402 0 : write(iulog,*) trim(errmsg)
403 0 : call endrun(errmsg)
404 : end if
405 : !
406 : ! n1_scha=19
407 : !
408 0 : istat = pio_inq_dimid(ncid, 'n1_scha', id)
409 0 : istat = pio_inquire_dimension(ncid, id, len=rd_n1_scha)
410 0 : if (rd_n1_scha /= n1_scha) then
411 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_n1_scha /= n1_scha rd_n1_scha = ', rd_n1_scha,' n1_scha = ', n1_scha
412 0 : write(iulog,*) trim(errmsg)
413 0 : call endrun(errmsg)
414 : end if
415 : !
416 : ! n2_scha=7
417 : !
418 0 : istat = pio_inq_dimid(ncid, 'n2_scha', id)
419 0 : istat = pio_inquire_dimension(ncid, id, len=rd_n2_scha)
420 0 : if (rd_n2_scha /= n2_scha) then
421 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_n2_scha /= n2_scha rd_n2_scha = ', rd_n2_scha,' n2_scha = ', n2_scha
422 0 : write(iulog,*) trim(errmsg)
423 0 : call endrun(errmsg)
424 : end if
425 : !
426 : ! n3_scha=68
427 : !
428 0 : istat = pio_inq_dimid(ncid, 'n3_scha', id)
429 0 : istat = pio_inquire_dimension(ncid, id, len=rd_n3_scha)
430 0 : if (rd_n3_scha /= n3_scha) then
431 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_n3_scha /= n3_scha rd_n3_scha = ', rd_n3_scha,' n3_scha = ', n3_scha
432 0 : write(iulog,*) trim(errmsg)
433 0 : call endrun(errmsg)
434 : end if
435 : !
436 : ! csize=28
437 : !
438 0 : istat = pio_inq_dimid(ncid, 'csize', id)
439 0 : istat = pio_inquire_dimension(ncid, id, len=rd_csize)
440 0 : if (rd_csize /= csize) then
441 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_csize /= csize rd_csize = ', rd_csize,' csize = ', csize
442 0 : write(iulog,*) trim(errmsg)
443 0 : call endrun(errmsg)
444 : end if
445 : !
446 : ! n_schfits=15
447 : !
448 0 : istat = pio_inq_dimid(ncid, 'n_schfits', id)
449 0 : istat = pio_inquire_dimension(ncid, id, len=rd_n_schfits)
450 0 : if (rd_n_schfits /= n_schfits) then
451 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_n_schfits /= n_schfits rd_n_schfits = ', &
452 0 : rd_n_schfits,' n_schfits = ', n_schfits
453 0 : write(iulog,*) trim(errmsg)
454 0 : call endrun(errmsg)
455 : end if
456 : !
457 : ! n_alschfits=18
458 : !
459 0 : istat = pio_inq_dimid(ncid, 'n_alschfits', id)
460 0 : istat = pio_inquire_dimension(ncid, id, len=rd_n_alschfits)
461 0 : if (rd_n_alschfits /= n_alschfits) then
462 0 : write(errmsg,"(a,i4,a,i4)") prefix//'rd_n_alschfits /= n_alschfits rd_n_alschfits = ',&
463 0 : rd_n_alschfits,' n_alschfits = ', n_alschfits
464 0 : write(iulog,*) trim(errmsg)
465 0 : call endrun(errmsg)
466 : end if
467 : !
468 : ! integer :: maxk_scha, maxm_scha, maxl_pot, maxm_pot
469 : ! maxk_scha = 18 ;
470 : ! maxm_scha = 6 ;
471 : ! maxl_pot = 12 ;
472 : ! maxm_pot = 2 ;
473 : !
474 0 : istat = pio_inq_dimid(ncid,"maxk_scha", id)
475 0 : istat = pio_inquire_dimension(ncid, id, len=maxk_scha)
476 0 : istat = pio_inq_dimid(ncid,"maxm_scha", id)
477 0 : istat = pio_inquire_dimension(ncid, id, len=maxm_scha)
478 0 : istat = pio_inq_dimid(ncid,"maxl_pot", id)
479 0 : istat = pio_inquire_dimension(ncid, id, len=maxl_pot)
480 0 : istat = pio_inq_dimid(ncid,"maxm_pot", id)
481 0 : istat = pio_inquire_dimension(ncid, id, len=maxm_pot)
482 :
483 : ! write(iulog,"('wei05sc: maxk_scha=',i3,' maxm_scha=',i3)") &
484 : ! maxk_scha,maxm_scha
485 : ! write(iulog,"('wei05sc: maxl_pot=',i3,' maxm_pot=',i3)") &
486 : ! maxl_pot,maxm_pot
487 : !
488 : ! Read variables:
489 : !
490 : ! double bndya(na):
491 0 : istat = pio_inq_varid(ncid, 'bndya', id)
492 0 : istat = pio_get_var(ncid, id,bndya)
493 : ! write(iulog,"('wei05sc: bndya=',/,(8f8.3))") bndya
494 : !
495 : ! double bndyb(nb):
496 0 : istat = pio_inq_varid(ncid, 'bndyb', id)
497 0 : istat = pio_get_var(ncid, id,bndyb)
498 : ! write(iulog,"('wei05sc: bndyb=',/,(8f8.3))") bndyb
499 : !
500 : ! double ex_bndy(nex):
501 0 : istat = pio_inq_varid(ncid, 'ex_bndy', id)
502 0 : istat = pio_get_var(ncid, id,ex_bndy)
503 : ! write(iulog,"('wei05sc: ex_bndy=',/,(8f8.3))") ex_bndy
504 : !
505 : ! double th0s(n3_scha):
506 0 : istat = pio_inq_varid(ncid, 'th0s', id)
507 0 : istat = pio_get_var(ncid, id,th0s)
508 : ! write(iulog,"('wei05sc: th0s=',/,(8f8.3))") th0s
509 : !
510 : ! double allnkm(n1_scha,n2_scha,n3_scha):
511 0 : istat = pio_inq_varid(ncid, 'allnkm', id)
512 0 : istat = pio_get_var(ncid, id,allnkm)
513 : ! write(iulog,"('wei05sc: allnkm min,max=',2e12.4)") minval(allnkm),maxval(allnkm)
514 : !
515 : ! int ab(csize):
516 0 : istat = pio_inq_varid(ncid, 'ab', id)
517 0 : istat = pio_get_var(ncid, id,ab)
518 : ! write(iulog,"('wei05sc: ab=',/,(10i4))") ab
519 : !
520 : ! int ls(csize):
521 0 : istat = pio_inq_varid(ncid, 'ls', id)
522 0 : istat = pio_get_var(ncid, id,ls)
523 : ! write(iulog,"('wei05sc: ls=',/,(10i4))") ls
524 : !
525 : ! int ms(csize):
526 0 : istat = pio_inq_varid(ncid, 'ms', id)
527 0 : istat = pio_get_var(ncid, id,ms)
528 : ! write(iulog,"('wei05sc: ms=',/,(10i4))") ms
529 : !
530 : ! double ex_epot(nex):
531 0 : istat = pio_inq_varid(ncid, 'ex_epot', id)
532 0 : istat = pio_get_var(ncid, id,ex_epot)
533 : ! write(iulog,"('wei05sc: ex_epot=',/,(8f8.3))") ex_epot
534 : !
535 : ! double ex_bpot(nex):
536 0 : istat = pio_inq_varid(ncid, 'ex_bpot', id)
537 0 : istat = pio_get_var(ncid, id,ex_bpot)
538 : ! write(iulog,"('wei05sc: ex_bpot=',/,(8f8.3))") ex_bpot
539 : !
540 : ! double epot_schfits(csize,n_schfits):
541 0 : istat = pio_inq_varid(ncid, 'epot_schfits', id)
542 0 : istat = pio_get_var(ncid, id,epot_schfits)
543 : ! write(iulog,"('wei05sc: epot_schfits min,max=',2e12.4)") &
544 : ! minval(epot_schfits),maxval(epot_schfits)
545 : !
546 : ! double bpot_schfits(csize,n_schfits):
547 0 : istat = pio_inq_varid(ncid, 'bpot_schfits', id)
548 0 : istat = pio_get_var(ncid, id,bpot_schfits)
549 : ! write(iulog,"('wei05sc: bpot_schfits min,max=',2e12.4)") &
550 : ! minval(bpot_schfits),maxval(bpot_schfits)
551 : !
552 : ! double epot_alschfits(csize,n_alschfits):
553 0 : istat = pio_inq_varid(ncid, 'epot_alschfits', id)
554 0 : istat = pio_get_var(ncid, id,epot_alschfits)
555 : ! write(iulog,"('wei05sc: epot_alschfits min,max=',2e12.4)") &
556 : ! minval(epot_alschfits),maxval(epot_alschfits)
557 : !
558 : ! double bpot_alschfits(csize,n_alschfits):
559 0 : istat = pio_inq_varid(ncid, 'bpot_alschfits', id)
560 0 : istat = pio_get_var(ncid, id,bpot_alschfits)
561 : ! write(iulog,"('wei05sc: bpot_alschfits min,max=',2e12.4)") &
562 : ! minval(bpot_alschfits),maxval(bpot_alschfits)
563 : !
564 : ! Close file:
565 0 : call cam_pio_closefile(ncid)
566 0 : if(masterproc) then
567 0 : write(iulog,"('wei05sc: completed read of file ',a)") trim(file)
568 : end if
569 :
570 0 : end subroutine read_wei05_ncfile
571 :
572 : !-----------------------------------------------------------------------
573 0 : subroutine setmodel(angle,bt,tilt,swvel,swden,model)
574 : !
575 : ! Calculate the complete set of the models' SCHA coeficients,
576 : ! given an aribitrary IMF angle (degrees from northward toward +Y),
577 : ! given byimf, bzimf, solar wind velocity (km/sec), and density.
578 : !
579 : ! Args:
580 : real(r8), intent(in) :: angle, bt, tilt, swvel, swden
581 : character(len=*), intent(in) :: model
582 : !
583 : ! Local:
584 : integer :: i, j
585 : real(r8) :: pi,stilt,stilt2,sw,swp,swe,c0,rang,cosa,sina,cos2a,sin2a
586 : real(r8) :: a(n_schfits)
587 : !
588 0 : if (trim(model) /= 'epot'.and.trim(model) /= 'bpot') then
589 0 : if (masterproc) then
590 0 : write(iulog, "('>>> model=',a)") trim(model)
591 : write(iulog, "(a)") &
592 0 : '>>> setmodel: model must be either ''epot'' or ''bpot'''
593 : end if
594 0 : call endrun("setmodel: model must be either 'epot' or 'bpot'")
595 : end if
596 : !
597 0 : pi = 4._r8 * atan(1._r8)
598 0 : rad2deg = 180._r8 / pi
599 0 : deg2rad = pi / 180._r8
600 : !
601 : ! write(iulog,"('setmodel call setboundary: model=',a,' swvel=',e12.4)") &
602 : ! model, swvel
603 :
604 0 : call setboundary(angle, bt, swvel, swden)
605 : !
606 0 : stilt = sin(tilt * deg2rad)
607 0 : stilt2 = stilt**2
608 0 : sw = bt * swvel/ 1000._r8
609 0 : if (trim(model) == 'epot') then
610 0 : swe = (1._r8-exp(-sw*ex_epot(2)))*sw**ex_epot(1)
611 : else
612 0 : swe = (1._r8-exp(-sw*ex_bpot(2)))*sw**ex_bpot(1)
613 : end if
614 0 : c0 = 1._r8
615 0 : swp = swvel**2 * swden*1.6726e-6_r8
616 0 : rang = angle*deg2rad
617 0 : cosa = cos(rang)
618 0 : sina = sin(rang)
619 0 : cos2a = cos(2._r8*rang)
620 0 : sin2a = sin(2._r8*rang)
621 0 : if (bt < 1._r8) then ! remove angle dependency for IMF under 1 nT
622 0 : cosa = -1._r8+bt*(cosa+1._r8)
623 0 : cos2a = 1._r8+bt*(cos2a-1._r8)
624 0 : sina = bt*sina
625 0 : sin2a = bt*sin2a
626 : end if
627 : a = (/c0, swe, stilt, stilt2, swp, &
628 : swe*cosa, stilt*cosa, stilt2*cosa, swp*cosa, &
629 : swe*sina, stilt*sina, stilt2*sina, swp*sina, &
630 0 : swe*cos2a, swe*sin2a/)
631 0 : if (trim(model) == 'epot') then
632 0 : esphc(:) = 0._r8
633 0 : do j=1,csize
634 0 : do i=1,n_schfits
635 0 : esphc(j) = esphc(j)+epot_schfits(i,j)*a(i)
636 : end do
637 : end do
638 : ! write(iulog,"('setmodel: esphc=',/,(6e12.4))") esphc
639 : else
640 0 : bsphc(:) = 0._r8
641 0 : do j=1,csize
642 0 : do i=1,n_schfits
643 0 : bsphc(j) = bsphc(j)+bpot_schfits(i,j)*a(i)
644 : end do
645 : end do
646 : ! write(iulog,"('setmodel: bsphc=',/,(6e12.4))") bsphc
647 : end if
648 0 : end subroutine setmodel
649 :
650 : !-----------------------------------------------------------------------
651 : !-----------------------------------------------------------------------
652 0 : subroutine wei05loc (ih, byimf, power, sunlon)
653 : ! ih=1,2 for SH,NH called from weimer05
654 : !
655 : ! (dimension 2 is for south, north hemispheres)
656 : ! Calculate offa, dskofa, rrad, phid, and phin from Weimer 2005 offc, dskofc, theta0
657 : ! Use Fig 8 of Heelis et al. [JGR, 85, 3315-3324, 1980]
658 : ! This shows: arad = 18.7 deg, crad = 16.7 deg (so arad = crad + 2 deg)
659 : ! offa = offc = 3 deg (so offa = offc)
660 : ! dskofc = 2 deg, dskofa = -0.5 deg (so dskofa = dskofc - 2.5 deg)
661 : ! Parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
662 : ! and phin (phin(MLT)=23.50 +/- 0.15By - 12)
663 : ! (In aurora_cons, phid=0., phin=180.*rtd)
664 : ! (For zero By, should be phid=21.39MLT*15*rtd, phin=11.5*15*rtd)
665 : ! 05/08: But formulae for ra-rc using IMF CP between 5-7 deg (not 2) so use
666 : ! difference of ra(max IMF CP or HP)-rc(IMF CP) as in aurora.F
667 : ! These are the dimensions and descriptions (corrected phid,n) from aurora.F:
668 : ! theta0(2), ! convection reversal boundary in radians
669 : ! offa(2), ! offset of oval towards 0 MLT relative to magnetic pole (rad)
670 : ! dskofa(2), ! offset of oval in radians towards 18 MLT (f(By))
671 : ! phid(2), ! dayside convection entrance in MLT-12 converted to radians (f(By))
672 : ! phid is the MLT-12 location of the cusp on the dayside
673 : ! phin(2), ! night convection entrance in MLT-12 converted to radians (f(By))
674 : ! rrad(2), ! radius of auroral circle in radians
675 : ! offc(2), ! offset of convection towards 0 MLT relative to mag pole (rad)
676 : ! dskofc(2) ! offset of convection in radians towards 18 MLT (f(By))
677 : ! sunlon: sun's longitude in dipole coordinates (see sub sunloc)
678 : !
679 : !
680 : ! Args:
681 : integer,intent(in) :: ih
682 : real(r8),intent(in) :: byimf
683 : real(r8),intent(in) :: power
684 : real(r8),intent(in) :: sunlon
685 : !
686 : ! Local:
687 : real(r8) :: rccp, racp, rahp, ramx, diffrac, plevel, tmltmin, tmltmax
688 : real(r8) :: offcdegp(2)
689 : integer :: i, j, j1, j2
690 : real(r8) :: vnx(2,2), hem, mltd, mltn
691 : integer :: inx(2,2)
692 : real(r8) :: offcdeg, dskof, arad, crad
693 : real(r8) :: byloc
694 :
695 : ! Limit size of byimf in phin and phid calculations (as in aurora.F)
696 : ! NOTE: This byloc is assymetric in hemisphere, which is probably not correct
697 0 : byloc = byimf
698 0 : if (byloc .gt. 7._r8) byloc = 7._r8
699 0 : if (byloc .lt. -11._r8) byloc = -11._r8
700 : !
701 : ! ih=1 is SH, ih=2 is NH
702 0 : if (ih .eq. 1) then
703 0 : j1 = 1
704 0 : j2 = nmlat0
705 0 : hem = -1._r8
706 : else
707 0 : j1 = nmlat0 + 1
708 0 : j2 = nmlat
709 0 : hem = 1._r8
710 : end if
711 : ! Print out un-revised values:
712 : ! write (6,"(1x,'Original convection/oval params (hem,By,off,dsk',
713 : ! | ',rad,phid,n=',10f9.4)") hem,byimf,offc(ih)*rtd,offa(ih)*rtd,
714 : ! | dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd,
715 : ! | phid(ih)*rtd/15.+12.,phin(ih)*rtd/15.+12.
716 : ! Find min/max
717 0 : vnx(ih,1) = 0._r8
718 0 : vnx(ih,2) = 0._r8
719 0 : do j=j1,j2
720 0 : do i=1,nmlonp1-1
721 0 : if (phihm(i,j) .gt. vnx(ih,2)) then
722 0 : vnx(ih,2) = phihm(i,j)
723 0 : inx(ih,2) = i
724 : end if
725 0 : if (phihm(i,j) .lt. vnx(ih,1)) then
726 0 : vnx(ih,1) = phihm(i,j)
727 0 : inx(ih,1) = i
728 : end if
729 : end do ! i=1,nmlonp1-1
730 : end do ! j=j1,j2
731 : ! 05/08: Calculate weictpoten in kV from Weimer model min/max in V
732 0 : weictpoten(ih) = 0.001_r8 * (vnx(ih,2) - vnx(ih,1))
733 0 : tmltmin = (ylonm(inx(ih,1))-sunlon) * r2d/15._r8 + 12._r8
734 : if (tmltmin > 24._r8) then
735 0 : tmltmin = tmltmin - 24._r8
736 : end if
737 0 : tmltmax = (ylonm(inx(ih,2))-sunlon) * r2d/15._r8 + 12._r8
738 : if (tmltmax > 24._r8) then
739 0 : tmltmax = tmltmax - 24._r8
740 : end if
741 : ! write (6,"('ih Bz By Hp ctpoten,wei min/max potV,lat,mlt=',i2,
742 : ! | 5f8.2,2x,e12.4,2f8.2,2x,e12.4,2f8.2))") ih,bzimf,byimf,power,
743 : ! | ctpoten,weictpoten(ih),
744 : ! | vnx(ih,1),ylatm(jnx(ih,1))*rtd,tmltmin,
745 : ! | vnx(ih,2),ylatm(jnx(ih,2))*rtd,tmltmax
746 : ! 05/08: From aurora_cons, calculate convection and aurora radii using IMF convection
747 : ! and power (plevel); racp (DMSP/NOAA) - rccp (AMIE) = 5.32 (Bz>0) to 6.62 (Bz<0) deg
748 : ! Heelis et al [1980, JGR, 85, pp 3315-3324] Fig 8: ra=rc+2deg, and is 2.5 deg to dusk
749 0 : rccp = -3.80_r8 + (8.48_r8*(weictpoten(ih)**0.1875_r8))
750 0 : racp = -0.43_r8 + (9.69_r8*(weictpoten(ih)**0.1875_r8))
751 0 : plevel = 0._r8
752 0 : if (power >= 1.00_r8) then
753 0 : plevel = 2.09_r8*log(power)
754 : end if
755 0 : rahp = 14.20_r8 + 0.96_r8*plevel
756 0 : ramx = max(racp, rahp)
757 0 : diffrac = ramx - rccp
758 :
759 : ! Set default values
760 : ! Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
761 : ! and phin (phin(MLT)=23.50 +/- 0.15By - 12)
762 0 : mltd = 9.39_r8 - hem*0.21_r8*byloc
763 0 : mltn = 23.50_r8 - hem*0.15_r8*byloc
764 0 : phid(ih) = (mltd-12._r8) * 15._r8 *d2r
765 0 : phin(ih) = (mltn-12._r8) * 15._r8 *d2r
766 : ! 05/18/08: Note that phid,phin are only for Heelis and are irrelevant for Weimer
767 : ! write (6,"(1x,'mltd mltn phid,n =',4f8.2)")
768 : ! | mltd,mltn,phid(ih)*rtd/15.,phin(ih)*rtd/15.
769 : ! Use default constant value of offcdegp from setboundary in Weimer 2005
770 0 : offcdeg = 4.2_r8
771 0 : offcdegp(ih) = offcdeg
772 0 : offc(ih) = offcdegp(ih) *d2r
773 0 : offa(ih) = offcdegp(ih) *d2r
774 : ! write (6,"(1x,'offcdeg,rad =',2e12.4)") offcdeg,offc(ih)
775 0 : dskof = 0._r8
776 0 : dskofc(ih) = dskof *d2r
777 : ! oval offset is 2.5 deg towards dawn (more neg dskof)
778 0 : dskofa(ih) = (dskof-2.5_r8) *d2r
779 : ! write (6,"(1x,'dskof,c,a=',3f8.2)")
780 : ! | dskof,dskofc(ih)*rtd,dskofa(ih)*rtd
781 : ! Set crad from bndyfitr/2 of setboundary of Weimer 2005
782 0 : crad = bndyfitr/2._r8
783 : ! write (6,"(1x,'wei05loc: ih,bz,y,crad =',i2,3f8.2)")
784 : ! | ih,bzimf,byimf,crad
785 : ! Fig 8 Heelis et al [1980]: ra=rc+2deg, and shifted 2.5 deg to dusk
786 0 : arad = crad + 2._r8
787 : ! 05/08: Make ra=rc+diffrac(=ramx-rccp) - same difference as in aurora.F
788 : ! Choose to have arad=crad(Weimer) + diffrac(same diff as in aurora.F)
789 0 : arad = crad + diffrac
790 : ! 08/08: OR make ra=ramx=max(racp,rahp) so diffrac=arad-crad
791 : ! diffrac2 = ramx - crad
792 : ! Choose to have arad=ramx (same as in aurora.F as determined by P/CP)
793 : ! arad = ramx
794 0 : theta0(ih) = crad *d2r
795 0 : rrad(ih) = arad *d2r
796 : ! write (6,"(1x,'radius: crad,rccp,racp,rahp diffa-c',
797 : ! | '(aurF,ramx-Weic) ramx,Weic+d,arad deg=',9f8.2)") crad,rccp,
798 : ! | racp,rahp,diffrac,diffrac2,ramx,crad+diffrac,arad
799 :
800 : ! Print out revised values (revised 05/08):
801 : ! write (6,"(1x,'Revised convection/oval params (off,dsk,',
802 : ! | 'rad,phid,n=',8f9.4)")offc(ih)*rtd,offa(ih)*rtd,
803 : ! | dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd,
804 : ! | phid(ih)*rtd/15.+12.,phin(ih)*rtd/15.+12.
805 :
806 0 : end subroutine wei05loc
807 :
808 : !-----------------------------------------------------------------------
809 : ! for now this is here ... might need to move to a gen util module
810 : !-----------------------------------------------------------------------
811 0 : function hp_from_bz_swvel(bz,swvel) result(hp)
812 : !
813 : ! Calculate hemispheric power from bz, swvel:
814 : ! Emery, et.al., (2008, in press, JGR)
815 : ! 6/3/08: Enforce minimum hp of 4.0 before *fac.
816 : ! 6/6/08: Reset minimum hp from 4.0 to 2.5 before *fac,
817 : ! as per Emery email of 6/5/08.
818 : !
819 : real(r8),intent(in) :: bz,swvel ! in
820 : real(r8) :: hp ! out
821 :
822 : real(r8), parameter :: fac = 2.0_r8
823 : !
824 0 : if (bz < 0._r8) then
825 0 : hp = 6.0_r8 + 3.3_r8*abs(bz) + (0.05_r8 + 0.003_r8*abs(bz))* (min(swvel,700._r8)-300._r8)
826 : else
827 0 : hp = 5.0_r8 + 0.05_r8 * (min(swvel,700._r8)-300._r8)
828 : end if
829 0 : hp = max(2.5_r8,hp)*fac
830 :
831 0 : end function hp_from_bz_swvel
832 :
833 : !-----------------------------------------------------------------------
834 0 : subroutine setboundary(angle,bt,swvel,swden)
835 : !
836 : ! Sets the coefficients that define the low-latitude boundary model,
837 : ! given the IMF and solar wind values.
838 : !
839 : implicit none
840 : !
841 : ! Args:
842 : real(r8),intent(in) :: angle,bt,swvel,swden
843 : !
844 : ! Local:
845 : integer :: i
846 : real(r8) :: swp,xc,theta,ct,st,cosa,btx,x(na),c(na)
847 : !
848 : ! Calculate the transformation matrix to the coordinate system
849 : ! of the offset pole.
850 : !
851 0 : xc = 4.2_r8
852 0 : theta = xc*(deg2rad)
853 0 : ct = cos(theta)
854 0 : st = sin(theta)
855 : !
856 0 : tmat(1,:) = (/ ct, 0._r8, st/)
857 0 : tmat(2,:) = (/ 0._r8, 1._r8, 0._r8/)
858 0 : tmat(3,:) = (/-st, 0._r8, ct/)
859 : !
860 : ! ttmat(1,:) = (/ct, 0._r8,-st/)
861 : ! ttmat(2,:) = (/ 0._r8,1._r8, 0._r8/)
862 : ! ttmat(3,:) = (/st, 0._r8, ct/)
863 : !
864 0 : swp = swden*swvel**2*1.6726e-6_r8 ! pressure
865 0 : cosa = cos(angle*deg2rad)
866 0 : btx = 1._r8-exp(-bt*ex_bndy(1))
867 0 : if (bt > 1._r8) then
868 0 : btx = btx*bt**ex_bndy(2)
869 : else
870 0 : cosa = 1._r8+bt*(cosa-1._r8) ! remove angle dependency for IMF under 1 nT
871 : end if
872 0 : x = (/1._r8, cosa, btx, btx*cosa, swvel, swp/)
873 0 : c = bndya
874 0 : bndyfitr = 0._r8
875 0 : do i=1,na
876 0 : bndyfitr = bndyfitr+x(i)*c(i)
877 :
878 : ! write(iulog,"('setboundry: i=',i3,' bndyfitr=',e12.4)") i,bndyfitr
879 :
880 : end do
881 0 : end subroutine setboundary
882 : !-----------------------------------------------------------------------
883 0 : subroutine epotval(lat,mlt,fill,epot)
884 : !
885 : ! Return the Potential (in kV) at given combination of def. latitude
886 : ! (lat) and MLT, in geomagnetic apex coordinates (practically identical
887 : ! to AACGM).
888 : ! If the location is outside of the model's low-latitude boundary, then
889 : ! the value "fill" is returned.
890 : !
891 : implicit none
892 : !
893 : ! Args:
894 : real(r8),intent(in) :: lat,mlt,fill
895 : real(r8),intent(out) :: epot
896 : !
897 : ! Local:
898 : integer :: inside,j,m,skip
899 : real(r8) :: z,phir,plm,colat,nlm
900 : real(r8) :: phim(2),cospm(2),sinpm(2)
901 : !
902 : ! checkinputs returns inside=1 if lat is inside model boundary,
903 : ! inside=0 otherwise. Phir and colat are also returned by checkinputs.
904 : !
905 0 : call checkinputs(lat,mlt,inside,phir,colat)
906 0 : if (inside == 0) then
907 0 : epot = fill
908 0 : return
909 : end if
910 : !
911 : ! IDL code:
912 : ! phim=phir # replicate(1,maxm) * ((indgen(maxm)+1) ## replicate(1,n_elements(phir)))
913 : ! where the '#' operator multiplies columns of first array by rows of second array,
914 : ! and the '##' operator multiplies rows of first array by columns of second array.
915 : ! Here, maxm == maxm_pot == 2, and phir is a scalar. The above IDL statement then
916 : ! becomes: phim = ([phir] # [1,1]) * ([1,2] ## [phir]) where phim will be
917 : ! dimensioned [1,2]
918 : !
919 0 : phim(1) = phir
920 0 : phim(2) = phir*2._r8
921 0 : cospm(:) = cos(phim(:))
922 0 : sinpm(:) = sin(phim(:))
923 : !
924 0 : z = 0._r8
925 0 : skip = 0 ! Added by B.Foster, 4/23/14
926 0 : do j=1,csize
927 0 : if (skip == 1) then
928 : skip = 0
929 : cycle
930 : end if
931 0 : m = ms(j)
932 0 : if (ab(j)==1) then
933 0 : plm = scplm(j,colat,nlm) ! scplm function is in this module
934 0 : skip = 0
935 0 : if (m == 0) then
936 0 : z = z+plm*esphc(j)
937 : else
938 0 : z = z+plm*(esphc(j)*cospm(m)+esphc(j+1)*sinpm(m))
939 0 : skip = 1
940 : end if
941 : end if ! ab(j)
942 : end do
943 0 : epot = z
944 : end subroutine epotval
945 : !-----------------------------------------------------------------------
946 0 : subroutine mpfac(lat,mlt,fill,fac)
947 : implicit none
948 : !
949 : ! Args:
950 : real(r8),intent(in) :: lat,mlt,fill
951 : real(r8),intent(out) :: fac
952 : !
953 : ! Local:
954 : integer :: j,m,inside,skip
955 : real(r8) :: phim(2),cospm(2),sinpm(2),cfactor
956 : real(r8) :: re,z,phir,plm,colat,nlm,pi
957 : !
958 0 : re = 6371.2_r8 + 110._r8 ! km radius (allow default ht=110)
959 : !
960 : ! checkinputs returns inside=1 if lat is inside model boundary,
961 : ! inside=0 otherwise. Phir and colat are also returned by checkinputs.
962 : !
963 0 : call checkinputs(lat,mlt,inside,phir,colat)
964 0 : if (inside == 0) then
965 0 : fac = fill
966 0 : return
967 : end if
968 : !
969 0 : phim(1) = phir
970 0 : phim(2) = phir*2._r8
971 0 : cospm(:) = cos(phim(:))
972 0 : sinpm(:) = sin(phim(:))
973 : !
974 0 : z = 0._r8
975 0 : skip = 0 ! Added by B.Foster, 4/23/14
976 0 : jloop: do j=1,csize
977 0 : if (skip == 1) then
978 : skip = 0
979 : cycle
980 : end if
981 0 : if (ls(j) >= 11) exit jloop
982 0 : m = ms(j)
983 0 : if (ab(j) == 1) then
984 0 : plm = scplm(j,colat,nlm) ! colat and nlm are returned (both reals)
985 0 : plm = plm*(nlm*(nlm+1._r8))
986 : !
987 : ! bsphc was calculated in setmodel (when setmodel called with 'bpot')
988 0 : if (m==0) then
989 0 : z = z-plm*bsphc(j)
990 : else
991 0 : z = z-(plm*(bsphc(j)*cospm(m)+bsphc(j+1)*sinpm(m)))
992 0 : skip = 1
993 : end if
994 : end if
995 : end do jloop ! j=1,csize
996 0 : pi = 4._r8*atan(1._r8)
997 0 : cfactor = -1.e5_r8/(4._r8*pi*re**2) ! convert to uA/m2
998 0 : z = z*cfactor
999 0 : fac = z
1000 : ! write(iulog,"('mpfac: lat=',f8.3,' mlt=',f8.3,' fac=',1pe12.4)") lat,mlt,fac
1001 : end subroutine mpfac
1002 : !-----------------------------------------------------------------------
1003 0 : real(r8) function scplm(index,colat,nlm)
1004 : !
1005 : ! Return Spherical Cap Harmonic Associated Legendre values, given colat
1006 : ! values and index i into array of L and M values.
1007 : !
1008 : implicit none
1009 : !
1010 : ! Args:
1011 : integer,intent(in) :: index
1012 : real(r8),intent(in) :: colat
1013 : real(r8),intent(out) :: nlm
1014 : !
1015 : ! Local:
1016 : integer :: i,j,l,m,skip
1017 : real(r8) :: th0,out(1),colata(1)
1018 : real(r8) :: cth(mxtablesize)
1019 : real(r8),save :: prevth0=1.e36_r8
1020 : integer,save :: tablesize
1021 : character(len=shr_kind_cl) :: errmsg
1022 : !
1023 0 : scplm = 0._r8
1024 0 : skip = 0 ! Added by B.Foster, 4/23/14
1025 0 : th0 = bndyfitr
1026 0 : if (prevth0 /= th0) then
1027 0 : tablesize = 3*nint(th0)
1028 0 : if (tablesize > mxtablesize) then
1029 : write(errmsg,"('>>> tablesize > mxtablesize: tablesize=',i8,' mxtablesize=',i8,' th0=',e12.4)") &
1030 0 : tablesize,mxtablesize,th0
1031 0 : write(iulog,*) trim(errmsg)
1032 0 : call endrun(errmsg)
1033 : end if
1034 0 : do i=1,tablesize
1035 0 : colattable(i) = real(i-1, r8) * (th0 / real(tablesize-1, r8))
1036 0 : cth(i) = cos(colattable(i) * deg2rad)
1037 : end do
1038 0 : prevth0 = th0
1039 0 : nlms = 0._r8 ! whole array init
1040 0 : do j=1,csize
1041 0 : if (skip == 1) then
1042 : skip = 0
1043 : cycle
1044 : end if
1045 0 : l = ls(j)
1046 0 : m = ms(j)
1047 0 : nlms(j) = nkmlookup(l,m,th0) ! nkmlookup in this module
1048 :
1049 : ! real :: plmtable(mxtablesize,csize)
1050 0 : call pm_n(m,nlms(j),cth,plmtable(1:tablesize,j),tablesize)
1051 0 : skip = 0
1052 0 : if (m /= 0 .and. ab(j) > 0) then
1053 0 : plmtable(1,j+1) = plmtable(1,j)
1054 0 : nlms(j+1) = nlms(j)
1055 0 : skip = 1
1056 : end if
1057 : end do ! j=1,csize
1058 : end if ! prevth0
1059 0 : nlm = nlms(index)
1060 0 : colata(1) = colat
1061 0 : call interpol_quad(plmtable(1:tablesize,index), &
1062 0 : colattable(1:tablesize),colata,out)
1063 0 : scplm = out(1)
1064 0 : end function scplm
1065 : !-----------------------------------------------------------------------
1066 0 : subroutine pm_n(m,r,cth,plmtable,tablesize)
1067 : !
1068 : ! Another SCHA function, returns the SCHA version of the associated
1069 : ! Legendre Polynomial, Pmn
1070 : !
1071 : ! Args:
1072 : integer,intent(in) :: m,tablesize
1073 : real(r8),intent(in) :: r
1074 : real(r8),intent(in) :: cth(tablesize)
1075 : real(r8),intent(out) :: plmtable(tablesize)
1076 : !
1077 : ! Local:
1078 : integer :: i,k
1079 : real(r8) :: rm,rk,div,ans,xn
1080 0 : real(r8),dimension(tablesize) :: a,x,tmp,table
1081 : !
1082 0 : if (m == 0) then
1083 0 : a = 1._r8 ! whole array op
1084 : else
1085 0 : do i=1,tablesize
1086 0 : a(i) = sqrt(1._r8-cth(i)**2)**m
1087 : end do
1088 : end if
1089 0 : xn = r*(r+1._r8)
1090 0 : x(:) = (1._r8-cth(:))/2._r8
1091 0 : table = a ! whole array init
1092 : k = 1
1093 : pmn_loop: do ! repeat-until loop in idl code
1094 0 : do i=1,tablesize
1095 0 : rm = real(m, r8)
1096 0 : rk = real(k, r8)
1097 0 : a(i) = a(i)*(x(i)*((rk+rm-1._r8)*(rk+rm)-xn)/(rk*(rk+rm)))
1098 0 : table(i) = table(i)+a(i) ! "result" in idl code
1099 : end do
1100 0 : k = k+1
1101 0 : do i=1,tablesize
1102 0 : div = abs(table(i))
1103 0 : if (div <= 1.e-6_r8) div = 1.e-6_r8
1104 0 : tmp(i) = abs(a(i)) / div
1105 : end do
1106 0 : if (maxval(tmp) < 1.e-6_r8) exit pmn_loop
1107 : end do pmn_loop
1108 0 : ans = km_n(m,r)
1109 :
1110 0 : plmtable(:) = table(:)*ans
1111 0 : end subroutine pm_n
1112 : !-----------------------------------------------------------------------
1113 0 : real(r8) function km_n(m,rn)
1114 : !
1115 : ! A normalization function used by the SCHA routines. See Haines.
1116 : !
1117 : ! Args:
1118 : integer,intent(in) :: m
1119 : real(r8),intent(in) :: rn
1120 : !
1121 : ! Local:
1122 : real(r8) :: rm
1123 : !
1124 0 : if (m == 0) then
1125 0 : km_n = 1._r8
1126 : return
1127 : end if
1128 0 : rm = real(m, r8)
1129 : km_n = sqrt(2._r8*exp(log_gamma(rn+rm+1._r8)-log_gamma(rn-rm+1._r8))) / &
1130 0 : (2._r8**m*factorial(m))
1131 0 : end function km_n
1132 : !-----------------------------------------------------------------------
1133 0 : real(r8) function nkmlookup(k, m, th0)
1134 : !
1135 : ! Given the size of a spherical cap, defined by the polar cap angle, th0,
1136 : ! and also the values of integers k and m, returns the value of n, a
1137 : ! real number (see Haines).
1138 : ! It uses interpolation from a lookup table that had been precomputed,
1139 : ! in order to reduce the computation time.
1140 : !
1141 : ! Args:
1142 : integer,intent(in) :: k,m
1143 : real(r8),intent(in) :: th0
1144 : !
1145 : ! Local:
1146 : integer :: kk,mm
1147 : real(r8) :: th0a(1),out(1)
1148 :
1149 0 : if (th0 == 90._r8) then
1150 0 : nkmlookup = real(k, r8)
1151 0 : return
1152 : end if
1153 0 : th0a(1) = th0
1154 0 : kk = k+1
1155 0 : mm = m+1
1156 0 : if (kk > maxk_scha) then
1157 0 : call interpol_quad(allnkm(maxk_scha,mm,:),th0s,th0a,out)
1158 : end if
1159 0 : if (mm > maxm_scha) then
1160 0 : call interpol_quad(allnkm(kk,maxm_scha,:),th0s,th0a,out)
1161 : end if
1162 0 : if (th0 < th0s(1)) then
1163 : write(iulog,"(a,e12.4,', th0s(1) = ',e12.4)") &
1164 0 : '>>> nkmlookup: th0 < th0s(1): th0 = ', th0, th0s(1)
1165 : end if
1166 0 : call interpol_quad(allnkm(kk,mm,:), th0s, th0a, out)
1167 0 : nkmlookup = out(1)
1168 0 : end function nkmlookup
1169 : !-----------------------------------------------------------------------
1170 0 : subroutine checkinputs(lat,mlt,inside,phir,colat)
1171 : !
1172 : ! Args:
1173 : real(r8), intent(in) :: lat,mlt
1174 : integer, intent(out) :: inside
1175 : real(r8), intent(out) :: phir,colat
1176 : !
1177 : ! Local:
1178 : real(r8) :: lon, tlat, tlon, radii
1179 : !
1180 0 : lon = mlt*15._r8
1181 0 : call dorotation(lat,lon,tlat,tlon)
1182 0 : radii = 90._r8-tlat
1183 0 : inside = 0
1184 0 : if (radii <= bndyfitr) then
1185 0 : inside = 1 ! bndyfitr from setboundary
1186 : end if
1187 0 : phir = tlon*deg2rad
1188 0 : colat = radii
1189 0 : end subroutine checkinputs
1190 : !-----------------------------------------------------------------------
1191 0 : subroutine dorotation(latin,lonin,latout,lonout)
1192 : !
1193 : ! Uses transformation matrices tmat and ttmat, to convert between
1194 : ! the given geomagnetic latatud/longitude, and the coordinate
1195 : ! system that is used within the model,that is offset from the pole.
1196 : !
1197 : ! Rotate Lat/Lon spherical coordinates with the transformation given
1198 : ! by saved matrix. The coordinates are assumed to be on a sphere of
1199 : ! Radius=1. Uses cartesian coordinates as an intermediate step.
1200 : !
1201 : implicit none
1202 : !
1203 : ! Args:
1204 : real(r8),intent(in) :: latin,lonin
1205 : real(r8),intent(out) :: latout,lonout
1206 : !
1207 : ! Local:
1208 : real(r8) :: latr,lonr,stc,ctc,sf,cf,a,b,pos(3)
1209 : integer :: i
1210 : !
1211 0 : latr = latin*deg2rad
1212 0 : lonr = lonin*deg2rad
1213 0 : stc = sin(latr)
1214 0 : ctc = cos(latr)
1215 0 : sf = sin(lonr)
1216 0 : cf = cos(lonr)
1217 0 : a = ctc*cf
1218 0 : b = ctc*sf
1219 : !
1220 : ! IDL code: Pos= TM ## [[A],[B],[STC]]
1221 : ! The ## operator multiplies rows of first array by columns of second array.
1222 : ! Currently, TM(3,3) = Tmat (or TTmat if "reversed" was set)
1223 : ! If called w/ single lat,lon, then a,b,stc are dimensioned (1), and
1224 : ! Pos is then (1,3)
1225 : !
1226 0 : do i=1,3
1227 0 : pos(i) = tmat(1,i)*a + tmat(2,i)*b + tmat(3,i)*stc
1228 : end do
1229 0 : latout = asin(pos(3))*rad2deg
1230 0 : lonout = atan2(pos(2),pos(1))*rad2deg
1231 0 : end subroutine dorotation
1232 : !-----------------------------------------------------------------------
1233 0 : subroutine interpol_quad(v,x,u,p)
1234 : !
1235 : ! f90 translation of IDL function interpol(v,x,u,/quadratic)
1236 : !
1237 : ! Args:
1238 : real(r8),intent(in) :: v(:),x(:),u(:)
1239 : real(r8),intent(out) :: p(:)
1240 : !
1241 : ! Local:
1242 : integer :: nv,nx,nu,i,ix
1243 : real(r8) :: x0,x1,x2
1244 : !
1245 0 : nv = size(v)
1246 0 : nx = size(x)
1247 0 : nu = size(u)
1248 0 : if (nx /= nv) then
1249 0 : p(:) = 0._r8
1250 : return
1251 : end if
1252 0 : do i = 1, nu
1253 0 : ix = value_locate(x,u(i))
1254 : ! 01/14 bae: interpol_quad in wei05sc.F is called when inside=1 or
1255 : ! radii<bndryfit for Weimer 2005. The fix to ix<=1 and ix>=nx
1256 : ! assures epot is non-zero near the pole (85.8mlat,0MLT) and
1257 : ! the boundary (bndryfit).
1258 0 : if (ix <=1) ix = 2
1259 0 : if (ix >=nx) ix = nx-1
1260 0 : x1 = x(ix)
1261 0 : x0 = x(ix-1)
1262 0 : x2 = x(ix+1)
1263 0 : p(i) = v(ix-1) * (u(i)-x1) * (u(i)-x2) / ((x0-x1) * (x0-x2)) + &
1264 0 : v(ix) * (u(i)-x0) * (u(i)-x2) / ((x1-x0) * (x1-x2)) + &
1265 0 : v(ix+1) * (u(i)-x0) * (u(i)-x1) / ((x2-x0) * (x2-x1))
1266 : end do
1267 : end subroutine interpol_quad
1268 : !-----------------------------------------------------------------------
1269 0 : integer function value_locate(vec,val)
1270 : !
1271 : ! f90 translation of IDL function value_locate
1272 : ! Return index i into vec for which vec(i) <= val >= vec(i+1)
1273 : ! Input vec must be monotonically increasing
1274 : !
1275 : implicit none
1276 : !
1277 : ! Args:
1278 : real(r8),intent(in) :: vec(:),val
1279 : !
1280 : ! Local:
1281 : integer :: n,i
1282 : !
1283 0 : value_locate = 0
1284 0 : n = size(vec)
1285 0 : if (val < vec(1)) then
1286 : return
1287 : end if
1288 0 : if (val > vec(n)) then
1289 0 : value_locate = n
1290 : return
1291 : end if
1292 0 : do i = 1, n-1
1293 0 : if (val >= vec(i) .and. val <= vec(i+1)) then
1294 0 : value_locate = i
1295 : return
1296 : end if
1297 : end do
1298 : end function value_locate
1299 : !-----------------------------------------------------------------------
1300 0 : real(r8) function factorial(n)
1301 : integer,intent(in) :: n
1302 : integer :: m
1303 0 : if (n <= 0) then
1304 0 : factorial = 0._r8
1305 : return
1306 : end if
1307 0 : if (n == 1) then
1308 0 : factorial = 1._r8
1309 : return
1310 : end if
1311 0 : factorial = real(n, r8)
1312 0 : do m = n-1,1,-1
1313 0 : factorial = factorial * real(m, r8)
1314 : end do
1315 : end function factorial
1316 : !-----------------------------------------------------------------------
1317 : !*********************** Copyright 1996,2001 Dan Weimer/MRC ***********************
1318 : ! COORDINATE TRANSFORMATION UTILITIES
1319 :
1320 : !NCAR Feb 01: Changed TRANS to GET_TILT s.t. the dipole tilt angle is
1321 : ! returned.
1322 :
1323 0 : real(r8) FUNCTION GET_TILT (YEAR,MONTH,DAY,HOUR)
1324 : ! SUBROUTINE TRANS(YEAR,MONTH,DAY,HOUR,IDBUG)
1325 : implicit none
1326 : real(r8) :: B3,B32,B33
1327 : integer :: IYR,JD,MJD,I,J,K
1328 : !NCAR
1329 :
1330 : INTEGER YEAR,MONTH,DAY,IDBUG
1331 : real(r8) :: HOUR
1332 : !
1333 : ! THIS SUBROUTINE DERIVES THE ROTATION MATRICES AM(I,J,K) FOR 11
1334 : ! TRANSFORMATIONS, IDENTIFIED BY K.
1335 : ! K=1 TRANSFORMS GSE to GEO
1336 : ! K=2 " GEO to MAG
1337 : ! K=3 " GSE to MAG
1338 : ! K=4 " GSE to GSM
1339 : ! K=5 " GEO to GSM
1340 : ! K=6 " GSM to MAG
1341 : ! K=7 " GSE to GEI
1342 : ! K=8 " GEI to GEO
1343 : ! K=9 " GSM to SM
1344 : ! K=10 " GEO to SM
1345 : ! K=11 " MAG to SM
1346 : !
1347 : ! IF IDBUG IS NOT 0, THEN OUTPUTS DIAGNOSTIC INFORMATION TO
1348 : ! FILE UNIT=IDBUG
1349 : !
1350 : INTEGER GSEGEO,GEOGSE,GEOMAG,MAGGEO
1351 : INTEGER GSEMAG,MAGGSE,GSEGSM,GSMGSE
1352 : INTEGER GEOGSM,GSMGEO,GSMMAG,MAGGSM
1353 : INTEGER GSEGEI,GEIGSE,GEIGEO,GEOGEI
1354 : INTEGER GSMSM,SMGSM,GEOSM,SMGEO,MAGSM,SMMAG
1355 :
1356 : PARAMETER (GSEGEO= 1,GEOGSE=-1,GEOMAG= 2,MAGGEO=-2)
1357 : PARAMETER (GSEMAG= 3,MAGGSE=-3,GSEGSM= 4,GSMGSE=-4)
1358 : PARAMETER (GEOGSM= 5,GSMGEO=-5,GSMMAG= 6,MAGGSM=-6)
1359 : PARAMETER (GSEGEI= 7,GEIGSE=-7,GEIGEO= 8,GEOGEI=-8)
1360 : PARAMETER (GSMSM = 9,SMGSM =-9,GEOSM =10,SMGEO=-10)
1361 : PARAMETER (MAGSM =11,SMMAG =-11)
1362 : !
1363 : ! The formal names of the coordinate systems are:
1364 : ! GSE - Geocentric Solar Ecliptic
1365 : ! GEO - Geographic
1366 : ! MAG - Geomagnetic
1367 : ! GSM - Geocentric Solar Magnetospheric
1368 : ! SM - Solar Magnetic
1369 : !
1370 : ! THE ARRAY CX(I) ENCODES VARIOUS ANGLES, STORED IN DEGREES
1371 : ! ST(I) AND CT(I) ARE SINES & COSINES.
1372 : !
1373 : ! Program author: D. R. Weimer
1374 : !
1375 : ! Some of this code has been copied from subroutines which had been
1376 : ! obtained from D. Stern, NASA/GSFC. Other formulas are from "Space
1377 : ! Physics Coordinate Transformations: A User Guide" by M. Hapgood (1991).
1378 : !
1379 : ! The formulas for the calculation of Greenwich mean sidereal time (GMST)
1380 : ! and the sun's location are from "Almanac for Computers 1990",
1381 : ! U.S. Naval Observatory.
1382 : !
1383 : real(r8) :: UT,T0,GMSTD,GMSTH,ECLIP,MA,LAMD,SUNLON
1384 :
1385 : !NCAR Feb 01: Eliminate unused routines from translib.for: ROTATE,
1386 : ! ROTATEV, FROMCART, TOCART, MLT, MAGLONG, SUNLOC. Remaining
1387 : ! are ADJUST and JULDAY
1388 : !NCAR Nov 02: Commons MFIELD and TRANSDAT now only in TRANS (GET_TILT)
1389 : ! So eliminate them as commons. For Fortran 90, eliminate
1390 : ! the DATA statement for assignments (not block_data)
1391 : ! COMMON/MFIELD/EPOCH,TH0,PH0,DIPOLE
1392 : ! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
1393 : !
1394 : real(r8) TH0,PH0 !,DIPOLE
1395 : real(r8) CX(9),ST(6),CT(6),AM(3,3,11)
1396 : !
1397 : ! TH0 = geog co-lat of NH magnetic pole
1398 : ! PH0 = geog longitude of NH magnetic pole
1399 : ! DIPOLE = magnitude of the B field in gauss at the equator
1400 :
1401 0 : TH0 = 11.19_r8
1402 0 : PH0 = -70.76_r8
1403 : ! DIPOLE = 0.30574_r8
1404 : !NCAR
1405 :
1406 : !NCAR Feb 01: Prevent debug printing to a file
1407 0 : IDBUG = 0
1408 : !NCAR
1409 :
1410 0 : IF(YEAR.LT.1900)THEN
1411 0 : IYR=1900+YEAR
1412 : ELSE
1413 0 : IYR=YEAR
1414 : END IF
1415 0 : UT=HOUR
1416 0 : JD=JULDAY(MONTH,DAY,IYR)
1417 0 : MJD=JD-2400001
1418 0 : T0=(real(MJD, r8) - 51544.5_r8) / 36525.0_r8
1419 : GMSTD=100.4606184_r8 + 36000.770_r8*T0 + 3.87933E-4_r8*T0*T0 + &
1420 0 : 15.0410686_r8*UT
1421 0 : CALL ADJUST(GMSTD)
1422 0 : GMSTH=GMSTD*24._r8/360._r8
1423 0 : ECLIP=23.439_r8 - 0.013_r8*T0
1424 0 : MA=357.528_r8 + 35999.050_r8*T0 + 0.041066678_r8*UT
1425 0 : CALL ADJUST(MA)
1426 0 : LAMD=280.460_r8 + 36000.772_r8*T0 + 0.041068642_r8*UT
1427 0 : CALL ADJUST(LAMD)
1428 0 : SUNLON=LAMD + (1.915_r8-0.0048_r8*T0)*SIND(MA) + 0.020_r8*SIND(2._r8*MA)
1429 0 : CALL ADJUST(SUNLON)
1430 : IF(IDBUG.NE.0)THEN
1431 : WRITE(IDBUG,*) YEAR,MONTH,DAY,HOUR
1432 : WRITE(IDBUG,*) 'MJD=',MJD
1433 : WRITE(IDBUG,*) 'T0=',T0
1434 : WRITE(IDBUG,*) 'GMSTH=',GMSTH
1435 : WRITE(IDBUG,*) 'ECLIPTIC OBLIQUITY=',ECLIP
1436 : WRITE(IDBUG,*) 'MEAN ANOMALY=',MA
1437 : WRITE(IDBUG,*) 'MEAN LONGITUDE=',LAMD
1438 : WRITE(IDBUG,*) 'TRUE LONGITUDE=',SUNLON
1439 : END IF
1440 :
1441 0 : CX(1)= GMSTD
1442 0 : CX(2) = ECLIP
1443 0 : CX(3) = SUNLON
1444 0 : CX(4) = TH0
1445 0 : CX(5) = PH0
1446 : ! Derived later:
1447 : ! CX(6) = Dipole tilt angle
1448 : ! CX(7) = Angle between sun and magnetic pole
1449 : ! CX(8) = Subsolar point latitude
1450 : ! CX(9) = Subsolar point longitude
1451 :
1452 0 : DO I=1,5
1453 0 : ST(I) = SIND(CX(I))
1454 0 : CT(I) = COSD(CX(I))
1455 : END DO
1456 : !
1457 0 : AM(1,1,GSEGEI) = CT(3)
1458 0 : AM(1,2,GSEGEI) = -ST(3)
1459 0 : AM(1,3,GSEGEI) = 0._r8
1460 0 : AM(2,1,GSEGEI) = ST(3)*CT(2)
1461 0 : AM(2,2,GSEGEI) = CT(3)*CT(2)
1462 0 : AM(2,3,GSEGEI) = -ST(2)
1463 0 : AM(3,1,GSEGEI) = ST(3)*ST(2)
1464 0 : AM(3,2,GSEGEI) = CT(3)*ST(2)
1465 0 : AM(3,3,GSEGEI) = CT(2)
1466 : !
1467 0 : AM(1,1,GEIGEO) = CT(1)
1468 0 : AM(1,2,GEIGEO) = ST(1)
1469 0 : AM(1,3,GEIGEO) = 0._r8
1470 0 : AM(2,1,GEIGEO) = -ST(1)
1471 0 : AM(2,2,GEIGEO) = CT(1)
1472 0 : AM(2,3,GEIGEO) = 0._r8
1473 0 : AM(3,1,GEIGEO) = 0._r8
1474 0 : AM(3,2,GEIGEO) = 0._r8
1475 0 : AM(3,3,GEIGEO) = 1._r8
1476 : !
1477 0 : DO I=1,3
1478 0 : DO J=1,3
1479 0 : AM(I,J,GSEGEO) = AM(I,1,GEIGEO)*AM(1,J,GSEGEI) + &
1480 0 : AM(I,2,GEIGEO)*AM(2,J,GSEGEI) + AM(I,3,GEIGEO)*AM(3,J,GSEGEI)
1481 : END DO
1482 : END DO
1483 : !
1484 0 : AM(1,1,GEOMAG) = CT(4)*CT(5)
1485 0 : AM(1,2,GEOMAG) = CT(4)*ST(5)
1486 0 : AM(1,3,GEOMAG) =-ST(4)
1487 0 : AM(2,1,GEOMAG) =-ST(5)
1488 0 : AM(2,2,GEOMAG) = CT(5)
1489 0 : AM(2,3,GEOMAG) = 0._r8
1490 0 : AM(3,1,GEOMAG) = ST(4)*CT(5)
1491 0 : AM(3,2,GEOMAG) = ST(4)*ST(5)
1492 0 : AM(3,3,GEOMAG) = CT(4)
1493 : !
1494 0 : DO I=1,3
1495 0 : DO J=1,3
1496 0 : AM(I,J,GSEMAG) = AM(I,1,GEOMAG)*AM(1,J,GSEGEO) + &
1497 0 : AM(I,2,GEOMAG)*AM(2,J,GSEGEO) + AM(I,3,GEOMAG)*AM(3,J,GSEGEO)
1498 : END DO
1499 : END DO
1500 : !
1501 0 : B32 = AM(3,2,GSEMAG)
1502 0 : B33 = AM(3,3,GSEMAG)
1503 0 : B3 = SQRT(B32*B32+B33*B33)
1504 0 : IF (B33.LE.0._r8) B3 = -B3
1505 : !
1506 0 : AM(2,2,GSEGSM) = B33/B3
1507 0 : AM(3,3,GSEGSM) = AM(2,2,GSEGSM)
1508 0 : AM(3,2,GSEGSM) = B32/B3
1509 0 : AM(2,3,GSEGSM) =-AM(3,2,GSEGSM)
1510 0 : AM(1,1,GSEGSM) = 1._r8
1511 0 : AM(1,2,GSEGSM) = 0._r8
1512 0 : AM(1,3,GSEGSM) = 0._r8
1513 0 : AM(2,1,GSEGSM) = 0._r8
1514 0 : AM(3,1,GSEGSM) = 0._r8
1515 : !
1516 0 : DO I=1,3
1517 0 : DO J=1,3
1518 0 : AM(I,J,GEOGSM) = AM(I,1,GSEGSM)*AM(J,1,GSEGEO) + &
1519 0 : AM(I,2,GSEGSM)*AM(J,2,GSEGEO) + AM(I,3,GSEGSM)*AM(J,3,GSEGEO)
1520 : END DO
1521 : END DO
1522 : !
1523 0 : DO I=1,3
1524 0 : DO J=1,3
1525 0 : AM(I,J,GSMMAG) = AM(I,1,GEOMAG)*AM(J,1,GEOGSM) + &
1526 0 : AM(I,2,GEOMAG)*AM(J,2,GEOGSM) + AM(I,3,GEOMAG)*AM(J,3,GEOGSM)
1527 : END DO
1528 : END DO
1529 : !
1530 0 : ST(6) = AM(3,1,GSEMAG)
1531 0 : CT(6) = SQRT(1._r8-ST(6)*ST(6))
1532 0 : CX(6) = ASIND(ST(6))
1533 :
1534 0 : AM(1,1,GSMSM) = CT(6)
1535 0 : AM(1,2,GSMSM) = 0._r8
1536 0 : AM(1,3,GSMSM) = -ST(6)
1537 0 : AM(2,1,GSMSM) = 0._r8
1538 0 : AM(2,2,GSMSM) = 1._r8
1539 0 : AM(2,3,GSMSM) = 0._r8
1540 0 : AM(3,1,GSMSM) = ST(6)
1541 0 : AM(3,2,GSMSM) = 0._r8
1542 0 : AM(3,3,GSMSM) = CT(6)
1543 : !
1544 0 : DO I=1,3
1545 0 : DO J=1,3
1546 0 : AM(I,J,GEOSM) = AM(I,1,GSMSM)*AM(1,J,GEOGSM) + &
1547 0 : AM(I,2,GSMSM)*AM(2,J,GEOGSM) + AM(I,3,GSMSM)*AM(3,J,GEOGSM)
1548 : END DO
1549 : END DO
1550 : !
1551 0 : DO I=1,3
1552 0 : DO J=1,3
1553 0 : AM(I,J,MAGSM) = AM(I,1,GSMSM)*AM(J,1,GSMMAG) + &
1554 0 : AM(I,2,GSMSM)*AM(J,2,GSMMAG) + AM(I,3,GSMSM)*AM(J,3,GSMMAG)
1555 : END DO
1556 : END DO
1557 : !
1558 0 : CX(7)=ATAN2D( AM(2,1,11) , AM(1,1,11) )
1559 0 : CX(8)=ASIND( AM(3,1,1) )
1560 0 : CX(9)=ATAN2D( AM(2,1,1) , AM(1,1,1) )
1561 :
1562 : IF(IDBUG.NE.0)THEN
1563 : WRITE(IDBUG,*) 'Dipole tilt angle=',CX(6)
1564 : WRITE(IDBUG,*) 'Angle between sun and magnetic pole=',CX(7)
1565 : WRITE(IDBUG,*) 'Subsolar point latitude=',CX(8)
1566 : WRITE(IDBUG,*) 'Subsolar point longitude=',CX(9)
1567 :
1568 : DO K=1,11
1569 : WRITE(IDBUG,1001) K
1570 : DO I=1,3
1571 : WRITE(IDBUG,1002) (AM(I,J,K),J=1,3)
1572 : END DO
1573 : END DO
1574 : 1001 FORMAT(' ROTATION MATRIX ',I2)
1575 : 1002 FORMAT(3F9.5)
1576 : END IF
1577 :
1578 : !NCAR Mar 96: return the dipole tilt from this function call.
1579 0 : GET_TILT = CX(6)
1580 : !NCAR
1581 :
1582 : RETURN
1583 : end function get_tilt
1584 : !-----------------------------------------------------------------------
1585 : !NCAR Feb 01: Eliminate unused routines from translib.for: ROTATE,
1586 : ! ROTATEV, FROMCART, TOCART, MLT, MAGLONG, SUNLOC. Remaining
1587 : ! are ADJUST and JULDAY
1588 : !NCAR
1589 0 : SUBROUTINE ADJUST(ANGLE)
1590 : implicit none
1591 : real(r8) :: angle
1592 : ! ADJUST AN ANGLE IN DEGREES TO BE IN RANGE OF 0 TO 360.
1593 : 10 CONTINUE
1594 0 : IF(ANGLE.LT.0._r8)THEN
1595 0 : ANGLE=ANGLE+360._r8
1596 0 : GOTO 10
1597 : END IF
1598 : 20 CONTINUE
1599 0 : IF(ANGLE.GE.360._r8)THEN
1600 0 : ANGLE=ANGLE-360._r8
1601 0 : GOTO 20
1602 : END IF
1603 0 : end subroutine adjust
1604 : !-----------------------------------------------------------------------
1605 0 : integer FUNCTION JULDAY(MM,ID,IYYY)
1606 : implicit none
1607 : integer :: igreg, iyyy, mm, id, jy, jm, ja
1608 : PARAMETER (IGREG=15+31*(10+12*1582))
1609 0 : IF (IYYY.EQ.0) call endrun('There is no Year Zero.')
1610 0 : IF (IYYY.LT.0) IYYY=IYYY+1
1611 0 : IF (MM.GT.2) THEN
1612 0 : JY=IYYY
1613 0 : JM=MM+1
1614 : ELSE
1615 0 : JY=IYYY-1
1616 0 : JM=MM+13
1617 : END IF
1618 0 : JULDAY=INT(365.25_r8*JY)+INT(30.6001_r8*JM)+ID+1720995
1619 0 : IF (ID+31*(MM+12*IYYY).GE.IGREG) THEN
1620 0 : JA=INT(0.01_r8*JY)
1621 0 : JULDAY=JULDAY+2-JA+INT(0.25_r8*JA)
1622 : END IF
1623 0 : end function julday
1624 : !-----------------------------------------------------------------------
1625 : SUBROUTINE CVT2MD(iulog,IYEAR,NDA,MON,DAY)
1626 : ! This sub converts NDA, the day number of the year, IYEAR,
1627 : ! into the appropriate month and day of month (integers)
1628 : implicit none
1629 : integer :: iulog,iyear,nda,mon,miss,numd,i
1630 : INTEGER DAY
1631 : INTEGER LMON(12)
1632 : PARAMETER (MISS=-32767)
1633 : SAVE LMON
1634 : DATA LMON/31,28,31,30,31,30,31,31,30,31,30,31/
1635 :
1636 : LMON(2)=28
1637 : IF(MOD(IYEAR,4) .EQ. 0)LMON(2)=29
1638 :
1639 : NUMD=0
1640 : DO 100 I=1,12
1641 : IF(NDA.GT.NUMD .AND. NDA.LE.NUMD+LMON(I))GO TO 200
1642 : NUMD=NUMD+LMON(I)
1643 : 100 CONTINUE
1644 : WRITE(iulog,'('' CVT2MD: Unable to convert year & day of year'', &
1645 : I5,'','',I5,''to month & day of month'')')IYEAR,NDA
1646 : MON = MISS
1647 : DAY = MISS
1648 : RETURN
1649 : 200 MON=I
1650 : DAY=NDA-NUMD
1651 : end subroutine cvt2md
1652 : !-----------------------------------------------------------------------
1653 : !
1654 : !NCAR Routines added to work around non-ANSI trig functions which
1655 : ! input degrees instead of radians: SIND, COSD, ASIND, ATAN2D
1656 :
1657 0 : FUNCTION SIND (DEG)
1658 : implicit none
1659 : real(r8) :: sind,d2r,r2d,deg
1660 : PARAMETER ( D2R = 0.0174532925199432957692369076847_r8 , &
1661 : R2D = 57.2957795130823208767981548147_r8)
1662 0 : SIND = SIN (DEG * D2R)
1663 0 : end function sind
1664 : !-----------------------------------------------------------------------
1665 0 : FUNCTION COSD (DEG)
1666 : implicit none
1667 : real(r8) :: cosd,d2r,r2d,deg
1668 : PARAMETER ( D2R = 0.0174532925199432957692369076847_r8 , &
1669 : R2D = 57.2957795130823208767981548147_r8)
1670 :
1671 0 : COSD = COS (DEG * D2R)
1672 0 : end function cosd
1673 : !-----------------------------------------------------------------------
1674 0 : FUNCTION ASIND (RNUM)
1675 : implicit none
1676 : real(r8) :: asind,d2r,r2d,rnum
1677 : PARAMETER ( D2R = 0.0174532925199432957692369076847_r8 , &
1678 : R2D = 57.2957795130823208767981548147_r8)
1679 0 : ASIND = R2D * ASIN (RNUM)
1680 0 : end function asind
1681 : !-----------------------------------------------------------------------
1682 0 : FUNCTION ATAN2D (RNUM1,RNUM2)
1683 : implicit none
1684 : real(r8) :: atan2d,d2r,r2d,rnum1,rnum2
1685 : PARAMETER ( D2R = 0.0174532925199432957692369076847_r8 , &
1686 : R2D = 57.2957795130823208767981548147_r8)
1687 0 : ATAN2D = R2D * ATAN2 (RNUM1,RNUM2)
1688 0 : end function atan2d
1689 : !-----------------------------------------------------------------------
1690 : end module wei05sc
|
