module efield
!---------------------------------------------------------------------
! description: calculates the electric potential for a given year,
! day of year,UT, F10.7, B_z(K_p)
! - low/midlatitudes electric potential is from an empirical model from
! L.Scherliess ludger@gaim.cass.usu.edu
! - high latitude electric potential is from Weimer96 model
! - the transition zone is smoothed
! - output is the horizontal global electric field in magnetic coordinates direction
! at every magnetic local time grid point expressed in degrees (0 deg-0MLT; 360deg 24 MLT)
!
! input
! integer :: iday, ! day number of year
! iyear ! year
! real:: ut, ! universal time
! F10.7, ! solar flux (see ionosphere module)
! bz ! component of IMF (see ionosphere module)
! output
! real :: &
! ed1(0:nmlon,0:nmlat), & ! zonal electric field Ed1 [V/m]
! ed2(0:nmlon,0:nmlat) ! meridional electric field Ed2/sin I_m [V/m]
!
! notes:
!
! - !to be done (commented out): input S_a F10.7/ Kp from WACCM and calculate B_z
! from these inputs
! - assume regular geomagnetic grid
! - uses average year 365.24 days/year 30.6001 day/mo s. Weimer
! - get_tilt works only for iyear >= 1900
! - Weimer model 1996, Dan Weimer (not with the updates from B.Emery)
! - fixed parameters: B_z, B_y units nT CHANGE THIS
! F10.7
! - we assume that the reference height is 300km for the emperical potential model
! - as a first approximation the electric field is constant in height
! WATCH what is the upper boundary condition in WACCM
! - for all the calculation done here we set the reference height to the same
! value as in tiegcm (hr=130km)
! - 12/15/03 input value iseasav : replaced by day -> month and day of month
! - 12/15/03 S_aM calculated according to Scherliess draft paper and added
! S_aM(corrected) = 90*(S_aM+1) to get variation in fig 1 Scherliess draft
!
! Apr 06 2012 Henry Juang, initial implement for nems
! Nov 20 2014 Jun Wang, change JULDAY to JULDAY_WAM
!
! Author: A. Maute Dec 2003 am 12/30/03
!------------------------------------------------------------------------------
! use shr_kind_mod, only: r8 => shr_kind_r8
! use physconst, only: pi
! use abortutils, only: endrun
! use cam_logfile, only: iulog
implicit none
public :: efield_init, ! interface routine
& get_efield ! interface routine
public :: ed1, ! zonal electric field Ed1 [V/m]
& ed2, ! meridional electric field Ed2 [V/m]
& potent, ! electric potential [V]
& nmlon, nmlat, ! dimension of mag. grid
& dlatm, dlonm, ! grid spacing of mag. grid
& ylonm, ylatm ! magnetic longitudes/latitudes (degc)
&,iday,iyear,iday_m,imo,f107d,by,bz,ut
public :: Coef , Cn,ML,MM1,MaxL,MaxM,MaxN,ALAMN,ALAMX,ALAMR,
&STPD,STP2,CSTP,SSTP,CX,ST,CT,AM,EPOCH,TH0,PH0,DIPOLE
! private
integer ::
& iday, ! day number of year
& iyear, ! year
& iday_m, ! day of month
& imo !month
real :: ut ! universal time
!----------------------------------------------------------------------
! solar parameters
!----------------------------------------------------------------------
real :: f107d ! 10.7 cm solar flux
real :: by ! By component of IMF [nT]
real :: bz ! Bz component of IMF [nT]
private
!----------------------------------------------------------------------
! mag. grid dimensions (assumed resolution of 2deg)
!----------------------------------------------------------------------
integer, parameter ::
&nmlon = 180, ! mlon
&nmlat = 90, ! mlat
&nmlath= nmlat/2, ! mlat/2
&nmlonh= nmlon/2, ! mlon/2
&nmlonp1 = nmlon+1, ! mlon+1
&nmlatp1 = nmlat+1, ! mlat+1
&iulog=10
real ::
& ylatm(0:nmlat), ! magnetic latitudes (deg)
& ylonm(0:nmlon), ! magnetic longitudes (deg)
& dlonm, ! delon lon grid spacing
& dlatm ! delat lat grid spacing
!----------------------------------------------------------------------
! array on magnetic grid:
!----------------------------------------------------------------------
real ::
& potent(0:nmlon,0:nmlat),! electric potential [V]
& ed1(0:nmlon,0:nmlat), ! zonal electric field Ed1 [V/m]
& ed2(0:nmlon,0:nmlat) ! meridional electric field Ed2/sin I_m [V/m]
real ::
& date, ! iyear+iday+ut
& day ! iday+ut
logical, parameter :: iutav=.false. ! .true. means UT-averaging
! .false. means no UT-averaging
! real, parameter :: v_sw = 400. ! solar wind velocity [km/s]
real, parameter :: v_sw = 450. ! solar wind velocity [km/s]
!----------------------------------------------------------------------
! boundary for Weimer
!----------------------------------------------------------------------
real, parameter :: bnd_wei = 44. ! colat. [deg]
integer :: nmlat_wei
!----------------------------------------------------------------------
! flag for choosing factors for empirical low latitude model
!----------------------------------------------------------------------
integer, parameter :: iseasav = 0 ! flag for season
!----------------------------------------------------------------------
! constants:
!----------------------------------------------------------------------
real, parameter ::
&r_e = 6.371e6, ! radius_earth [m] (same as for apex.F90)
& h_r = 130.0e3, ! reference height [m] (same as for apex.F90)
&dy2yr= 365.24, ! day per avg. year used in Weimer
&dy2mo= 30.6001, ! day per avg. month used in Weimer
&pi=3.141592653
real
& rtd , ! radians -> deg
& dtr, ! deg -> radians
& sqr2,
& hr2rd, ! pi/12 hrs
& dy2rd, ! 2*pi/365.24 average year
& deg2mlt, ! for mlon to deg
& mlt2deg, ! for mlt to mlon
& sinIm_mag(0:nmlat) ! sinIm
integer :: jmin, jmax ! latitude index for interpolation of
! northward e-field ed2 at mag. equator
!----------------------------------------------------------------------
! for spherical harmonics
!----------------------------------------------------------------------
integer, parameter ::
& nm = 19,
& mm = 18,
& nmp = nm + 1,
& mmp = mm + 1
real :: r(0:nm,0:mm) ! R_n^m
real :: pmopmmo(0:mm) ! sqrt(1+1/2m)
!----------------------------------------------------------------------
! index for factors f_m(mlt),f_l(UT),f_-k(d)
!----------------------------------------------------------------------
integer, parameter :: ni = 1091 ! for n=12 m=-18:18
integer :: imax ! max number of index
integer,dimension(0:ni) :: kf,lf, mf, nf, jf
real :: ft(1:3,0:2) ! used for f_-k(season,k)
real :: a_klnm(0:ni) ! A_klm
real :: a_lf(0:ni) ! A_klmn^lf for minimum
real :: a_hf(0:ni) ! A_klmn^hf for maximum
!----------------------------------------------------------------------
!replace wei96.f common block
!----------------------------------------------------------------------
real :: Coef(0:1,0:8,0:3)
real :: Cn(0:3,0:1,0:4,0:1,0:8,0:3)
integer ML,MM1,MaxL,MaxM,MaxN
real ALAMN,ALAMX,ALAMR,STPD,STP2,CSTP,SSTP
real CX(9),ST(6),CT(6),AM(3,3,11)
real , parameter :: EPOCH=1980.,TH0=11.19,PH0=-70.76,
& DIPOLE=.30574
!----------------------------------------------------------------------
! high_latitude boundary
!----------------------------------------------------------------------
real, parameter ::
&ef_max = 0.015, ! max e-field for high latitude boundary location [V/m]
&lat_sft = 54. ! shift of highlat_bnd to 54 deg
integer :: ilat_sft ! index of shift for high latitude boundary
integer, parameter :: nmax_sin = 2 ! max. wave number to be represented
logical, parameter :: debug =.false.
!
contains
subroutine efield_init
!hmhj subroutine efield_init(efield_lflux_file, efield_hflux_file,
!hmhj&efield_wei96_file)
!--------------------------------------------------------------------
! Purpose: read in and set up coefficients needed for electric field
! calculation (independent of time & geog. location)
!
! Method:
!
! Author: A. Maute Dec 2003 am 12/17/03
!-------------------------------------------------------------------
!hmhj character(len=*), intent(in) :: efield_lflux_file
!hmhj character(len=*), intent(in) :: efield_hflux_file
!hmhj character(len=*), intent(in) :: efield_wei96_file
character(len=*), parameter ::
& efield_lflux_file='global_idea_coeff_lflux.dat',
& efield_hflux_file='global_idea_coeff_hflux.dat',
& efield_wei96_file='global_idea_wei96.cofcnts'
call constants ! calculate constants
!-----------------------------------------------------------------------
! low/midlatitude potential from Scherliess model
!-----------------------------------------------------------------------
call read_acoef (efield_lflux_file, efield_hflux_file) ! read in A_klnm for given S_aM
call index_quiet ! set up index for f_m(mlt),f_l(UT),f_-k(d)
call prep_fk ! set up the constant factors for f_k
call prep_pnm ! set up the constant factors for P_n^m & dP/d phi
!-----------------------------------------------------------------------
!following part should be independent of time & location if IMF constant
!-----------------------------------------------------------------------
call ReadCoef (efield_wei96_file)
end subroutine efield_init
subroutine get_efield
!-----------------------------------------------------------------------
! Purpose: calculates the global electric potential field on the
! geomagnetic grid (MLT in deg) and derives the electric field
!
! Method:
!
! Author: A. Maute Dec 2003 am 12/17/03
!-----------------------------------------------------------------------
! use time_manager, only : get_curr_calday, get_curr_date
! use mo_solar_parms, only : get_solar_parms
! use mag_parms, only : get_mag_parms
! use cam_control_mod, only: magfield_fix_year
! use spmd_utils, only: masterproc
integer :: idum1, idum2, tod ! time of day [s]
real kp
!-----------------------------------------------------------------------
! get current calendar day of year & date components
! valid at end of current timestep
!-----------------------------------------------------------------------
! iday = get_curr_calday() ! day of year
! call get_curr_date (iyear,imo,iday_m,tod)! year, time of day [sec]
! iyear = magfield_fix_year
! iyear = 1995
! if( iyear < 1900 ) then
! write(iulog,"(/,'>>> get_efield: year < 1900 not possible:
! &year=',i5)") iyear
! call endrun
! end if
tod=ut*3600.
! ut = tod/3600. ! UT of day [sec]
!-----------------------------------------------------------------------
! get solar parms
!-----------------------------------------------------------------------
! call get_solar_parms( f107_s = f107d )
!-----------------------------------------------------------------------
! get mag parms
!-----------------------------------------------------------------------
! call get_mag_parms( by = by, bz = bz )
! print*,by,bz,f107d,ut
!#ifdef EFIELD_DIAGS
! if( masterproc ) then
! write(iulog,*) 'get_efield: f107d,by,bz = ', f107d,by,bz
! end if
!#endif
!-----------------------------------------------------------------------
! ajust S_a
!-----------------------------------------------------------------------
call adj_S_a
!-----------------------------------------------------------------------
! calculate global electric potential
!-----------------------------------------------------------------------
call GlobalElPotential
! print*,'pot_efield',potent(149,66),potent(149,64)
!-----------------------------------------------------------------------
! calculate derivative of global electric potential
!-----------------------------------------------------------------------
call DerivPotential
! print*,'ed2_efield',ed2(149,65),potent(149,66),potent(149,64)
end subroutine get_efield
subroutine GlobalElPotential
!-----------------------------------------------------------------------
! Purpose: calculates the global electric potential field on the
! geomagnetic grid (MLT in deg)
!
! Method: rewritten code from Luedger Scherliess (11/20/99 LS)
! routine to calculate the global electric potential in magnetic
! Apex coordinates (Latitude and MLT).
! High Latitude Model is Weimer 1996.
! Midlatitude model is Scherliess 1999.
! Interpolation in a transition region at about 60 degree
! magnetic apex lat
!
! Author: A. Maute Dec 2003 am 12/17/03
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer :: ilon, ilat, idlat
integer :: ihlat_bnd(0:nmlon) ! high latitude boundary
integer :: itrans_width(0:nmlon) ! width of transition zone
real :: mlt, mlon, mlat, mlat_90, pot
real :: pot_midlat(0:nmlon,0:nmlat) ! potential from L. Scherliess model
real :: pot_highlat(0:nmlon,0:nmlat) ! potential from Weimer model
real :: pot_highlats(0:nmlon,0:nmlat)! smoothed potential from Weimer model
!-----------------------------------------------------------------------
! Externals
!-----------------------------------------------------------------------
real,external :: EpotVal ! in wei96.f
!-----------------------------------------------------------------------
! convert to date and day
!-----------------------------------------------------------------------
day = iday + ut/24.
date = iyear + day/dy2yr
!-----------------------------------------------------------------------
! low/midlatitude electric potential - empirical model Scherliess 1999
!-----------------------------------------------------------------------
!$omp parallel do private(ilat, ilon, mlat, pot)
do ilat = 0,nmlath ! Calculate only for one magn. hemisphere
mlat = ylatm(ilat) ! mag. latitude
do ilon = 0,nmlon ! lon. loop
call efield_mid( mlat, ylonm(ilon), pot )
pot_midlat(ilon,ilat+nmlath) = pot ! SH/NH symmetry
pot_midlat(ilon,nmlath-ilat) = pot
end do
end do
! print*,'www1','midlat',pot_midlat(149,66)
!-----------------------------------------------------------------------
! hight latitude potential from Weimer model
! at the poles Weimer potential is not longitudinal dependent
!-----------------------------------------------------------------------
call prep_weimer ! calculate IMF angle & magnitude, tilt
!$omp parallel do private(ilat, ilon, mlat_90, pot)
do ilat = 0,nmlat_wei ! Calculate only for one magn. hemisphere
mlat_90 = 90. - ylatm(ilat) ! mag. latitude
do ilon = 0,nmlon
pot = 1000.*EpotVal( mlat_90, ylonm(ilon)*deg2mlt ) ! calculate potential (kv -> v)
!-----------------------------------------------------------------------
! NH/SH symmetry
!-----------------------------------------------------------------------
pot_highlat(ilon,ilat) = pot
pot_highlat(ilon,nmlat-ilat) = pot
pot_highlats(ilon,ilat) = pot
pot_highlats(ilon,nmlat-ilat) = pot
! bad value com from EpotVal
! if(ilat.eq.22.and.ilon.eq.148)
! & print*,'www2',ilat,ilon,pot,mlat_90,ylonm(ilon)*deg2mlt
end do
end do
! print*,'www2','highlat',ut,by,bz,pot_highlat(0:180,68),nmlat_wei
!-----------------------------------------------------------------------
! weighted smoothing of high latitude potential
!-----------------------------------------------------------------------
idlat = 2 ! smooth over -2:2 = 5 grid points
call pot_latsmo( pot_highlats, idlat )
! print*,'www2','highlat',ut,pot_highlat(0:180,45)
!-----------------------------------------------------------------------
! calculate the height latitude bounday ihl_bnd
! 1. calculate E field from weimar model
! boundary is set where the total electric field exceeds
! 0.015 V/m (corresp. approx. to 300 m/s)
! 2. moved halfways to 54 deg
! output : index 0-pole nmlath-equator
!-----------------------------------------------------------------------
call highlat_getbnd( ihlat_bnd )
!-----------------------------------------------------------------------
! 3. adjust high latitude boundary sinusoidally
! calculate width of transition zone
!-----------------------------------------------------------------------
call bnd_sinus( ihlat_bnd, itrans_width )
!-----------------------------------------------------------------------
! 4. ajust high latitude potential to low latitude potential
!-----------------------------------------------------------------------
! print*,'www30',ihlat_bnd
call highlat_adjust( pot_highlats, pot_highlat, pot_midlat,
&ihlat_bnd )
! print*,'www3','highlat',ut,pot_highlat(145:153,68)
! print*,'www3','midlat',ut,pot_midlat(145:153,68)
!-----------------------------------------------------------------------
! interpolation of high and low/midlatitude potential in the
! transition zone and put it into global potent array
!-----------------------------------------------------------------------
call interp_poten( pot_highlats, pot_highlat, pot_midlat,
&ihlat_bnd, itrans_width)
! print*,'www4','potent',ut,by,bz,potent(0:181,68)
!-----------------------------------------------------------------------
! potential weighted smoothing in latitude
!-----------------------------------------------------------------------
idlat = 2 ! smooth over -2:2 = 5 grid points
call pot_latsmo2( potent, idlat )
! print*,'www5','pot_efield',potent(149,68)
!-----------------------------------------------------------------------
! potential smoothing in longitude
!-----------------------------------------------------------------------
idlat = nmlon/48 ! smooth over -idlat:idlat grid points
call pot_lonsmo( potent, idlat )
! print*,'www6','pot_efield',ut,by,bz,potent(0:180,68)
!-----------------------------------------------------------------------
! output
!-----------------------------------------------------------------------
! output ( change later to netcdf file)
! do ilat=0,nmlat
! do ilon=0,nmlon
! write(iulog,'(4(x,f12.5))') ylatm(ilat),ylonm(ilon), &
! potent(ilon,ilat),potent(ilon,nmlat-ilat)
! write(iulog,'(4(x,f12.5))') ylatm(ilat),ylonm(ilon), &
! potent(ilon,ilat),potent(ilon,nmlat-ilat)
! write(iulog,'(f10.3)') potent(ilon,ilat)
! end do
! end do
end subroutine GlobalElPotential
subroutine ff( ph, mt, f )
!-----------------------------------------------------------------------
!Purpose: calculate F for normalized associated Legendre polynomial P_n^m
! Ref.: Richmond J.Atm.Ter.Phys. 1974
!
! Method: f_m(phi) = sqrt(2) sin(m phi) m > 0
! = 1 m = 0
! = sqrt(2) cos(m phi) m < 0
!
! Author: A. Maute Nov 2003 am 11/18/03
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
! dummy arguments
!-----------------------------------------------------------------------
integer,intent(in) :: mt
real,intent(in) :: ph ! geo. longitude of 0SLT (ut*15)
real,intent(out) :: f(-mt:mt)
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer :: m, i, j, mmo
real :: sp, cp
sp = sin( ph/rtd )
cp = cos( ph/rtd )
f(0) = 1.e0
f(-1) = sqr2*cp
f(1) = sqr2*sp
do m = 2,mt
mmo = m - 1
f(m) = f(-mmo)*sp + cp*f(mmo)
f(-m) = f(-mmo)*cp - sp*f(mmo)
end do
end subroutine ff
subroutine pnm( ct, p )
!----------------------------------------------------------------
! Purpose: normalized associated Legendre polynomial P_n^m
! Ref.: Richmond J.Atm.Ter.Phys. 1974
! Method:
! P_m^m = sqrt(1+1/2m)*si*P_m-1^m-1 m>0
! P_n^m = [cos*P_n-1^m - R_n-1^m*P_n-2^m ]/R_n^m n>m>=0
! dP/d phi = n*cos*P_n^m/sin-(2*n+1)*R_n^m*P_n-1^m/sin n>=m>=0
! R_n^m = sqrt[ (n^2-m^2)/(4n^2-1) ]
!
! Author: A. Maute Nov 2003 am 11/18/03
!--------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
! dummy arguments
!-----------------------------------------------------------------------
real, intent(inout) :: ct ! cos(colat)
real, intent(inout) :: p(0:nm,0:mm)
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer :: mp, m, n, np
real :: pm2, st
! ct = min( ct,.99 ) ! cos(colat)
st = sqrt( 1. - ct*ct ) ! sin(colat)
p(0,0) = 1.
do mp = 1,mmp ! m+1=1,mm+1
m = mp - 1
if( m >= 1 ) then
p(m,m) = pmopmmo(m)*p(m-1,m-1)*st
end if
pm2 = 0.
do n = mp,nm ! n=m+1,N
np = n + 1
p(n,m) = (ct*p(n-1,m) - r(n-1,m)*pm2)/r(n,m)
pm2 = p(n-1,m)
end do
end do
end subroutine pnm
subroutine prep_pnm
!-----------------------------------------------------------------
! Purpose: constant factors for normalized associated Legendre polynomial P_n^m
! Ref.: Richmond J.Atm.Ter.Phys. 1974
!
! Method:
! PmoPmmo(m) = sqrt(1+1/2m)
! R_n^m = sqrt[ (n^2-m^2)/(4n^2-1) ]
!
! Author: A. Maute Nov 2003 am 11/18/03
!-----------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer :: mp, m, n
real :: xms, xns, den
do mp = 1, mmp ! m+1 = 1,mm+1
m = mp - 1
xms = m*m
if( mp /= 1 ) then
pmopmmo(m) = sqrt( 1. + .5/M )
end if
do n = m,nm ! n = m,N
xns = n*n
den = max(4.*xns - 1.,1.)
r(n,m) = sqrt( (xns - xms)/den )
end do
end do
end subroutine prep_pnm
subroutine index_quiet
!-----------------------------------------------------------------
! Purpose: set up index for factors f_m(mlt),f_l(UT),f_-k(d) to
! describe the electric potential Phi for the empirical model
!
! Method:
! Phi = sum_k sum_l sum_m sum_n [ A_klmn * P_n^m *f_m(mlt)*f_l(UT)*f_-k(d)]
! - since the electric potential is symmetric about the equator
! n+m odd terms are set zero resp. not used
! - in the summation for calculation Phi the index have the following
! range n=1,12 and m=-n,n, k=0,2 l=-2,2
!
! Author: A. Maute Nov 2003 am 11/18/03
!----------------------------------------------------------------
implicit none
!----------------------------------------------------------------
! ... local variables
!----------------------------------------------------------------
integer :: i, j, k, l, n, m
i = 0 ! initialize
j = 1
do k = 2,0,-1
do l = -2,2
if( k == 2 .and. abs(l) == 2 ) then
cycle
end if
do n = 1,12
do m = -18,18
if( abs(m) <= n ) then ! |m| < n
if( (((n-m)/2)*2) == (n-m) ) then ! only n+m even
if( n-abs(m) <= 9 ) then ! n-|m| <= 9 why?
kf(i) = 2-k
lf(i) = l
nf(i) = n
mf(i) = m
jf(i) = j
i = i + 1 ! counter
end if
end if
end if
end do ! m
end do ! n
end do ! l
end do ! k
imax = i - 1
if(imax /= ni ) then ! check if imax == ni
! write(iulog,'(a19,i5,a18,i5)') 'index_quiet: imax= ',imax,
! & ' not equal to ni =',ni
stop
end if
! if(debug) write(iulog,*) 'imax=',imax
end subroutine index_quiet
subroutine read_acoef (efield_lflux_file, efield_hflux_file)
!----------------------------------------------------------------
! Purpose:
! 1. read in coefficients A_klmn^lf for solar cycle minimum and
! A_klmn^hf for maximum
! 2. adjust S_a (f107d) such that if S_a<80 or S_a > 220 it has reasonable numbers
! S_aM = [atan{(S_a-65)^2/90^2}-a90]/[a180-a90]
! a90 = atan [(90-65)/90]^2
! a180 = atan [(180-65)/90]^2
! 3. inter/extrapolation of the coefficient to the actual flux which is
! given by the user
! A_klmn = S_aM [A_klmn^hf-A_klmn^lf]/90. + 2*A_klmn^lf-A_klmn^hf
!
! Method:
!
! Author: A. Maute Nov 2003 am 11/19/03
!---------------------------------------------------------------
! use ioFileMod, only : getfil
! use units, only : getunit, freeunit
character(len=*), intent(in) :: efield_lflux_file
character(len=*), intent(in) :: efield_hflux_file
integer :: i,ios,unit,istat
! character (len=13):: locfn
character (len=256):: locfn
!------------------------------------------------------------------
! get coefficients file for solar minimum:
!-----------------------------------------------------------------
! unit = getunit()
unit = 11
! call getfil( efield_lflux_file, locfn, 0 )
locfn=efield_lflux_file
!------------------------------------------------------------------
! open datafile with coefficients A_klnm
!------------------------------------------------------------------
! write(iulog,*) 'read_acoef: open file ',trim(locfn),
! &' unit ',unit
open(unit=unit,file=trim(locfn),
& status = 'old',iostat = ios)
! if(ios.gt.0) then
! write(iulog,*)
! &'read_acoef: error in opening coeff_lf file',
! &' unit ',unit
! call endrun
! end if
!----------------------------------------------------------------------------
! read datafile with coefficients A_klnm
!--------------------------------------------------------------------
! write(iulog,*) 'read_acoef: read file ',trim(locfn),' unit ',unit
read(unit,*,iostat = ios) a_lf
! if(ios.gt.0) then
! write(iulog,*)
! &'read_acoef: error in reading coeff_lf file',' unit ',unit
! call endrun
! end if
!--------------------------------------------------------------------
! close & free unit
!--------------------------------------------------------------------
close(unit)
! call freeunit(unit)
! write(iulog,*) 'read_acoef: free unit ',unit
!--------------------------------------------------------------------
! get coefficients file for solar maximum:
!--------------------------------------------------------------------
! unit = getunit()
unit = 10
! call getfil( efield_hflux_file, locfn, 0 )
locfn= efield_hflux_file
!-------------------------------------------------------------------
! open datafile with coefficients A_klnm
!------------------------------------------------------------------
! write(iulog,*) 'read_acoef: open file ',trim(locfn),' unit ',unit
open(unit=unit,file=trim(locfn),
& status = 'old',iostat = ios)
! if(ios.gt.0) then
! write(iulog,*)
! &'read_acoef: error in opening coeff_hf file',' unit ',unit
! call endrun
! end if
!-----------------------------------------------------------------
! read datafile with coefficients A_klnm
!----------------------------------------------------------------
! write(iulog,*) 'read_acoef: read file ',trim(locfn)
read(unit,*,iostat = ios) a_hf
! if(ios.gt.0) then
! write(iulog,*)
! &'read_acoef: error in reading coeff_hf file',' unit ',unit
! call endrun
! end if
!---------------------------------------------------------------
! close & free unit
!--------------------------------------------------------------
close(unit)
! call freeunit(unit)
! write(iulog,*) 'read_acoef: free unit ',unit
end subroutine read_acoef
subroutine adj_S_a
!------------------------------------------------------------------
! adjust S_a -> S_aM eqn.8-11 Scherliess draft
!------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------
! local variables
!------------------------------------------------------------------
integer :: i
real :: x2, y2, a90, a180, S_aM
x2 = 90.*90.
y2 = (90. - 65.)
y2 = y2*y2
a90 = atan2(y2,x2)
y2 = (180. - 65.)
y2 = y2*y2
a180 = atan2(y2,x2)
! y2 = (S_a-65.)
y2 = (f107d - 65.)
y2 = y2*y2
S_aM = (atan2(y2,x2) - a90)/(a180 - a90)
S_aM = 90.*(1. + S_aM)
! if(debug) write(iulog,*) 'f107d=',f107d,' S_aM =',S_aM
! if(debug) write(iulog,*) 'By=',by
!-----------------------------------------------------------------
! inter/extrapolate to S_a (f107d)
!----------------------------------------------------------------
do i = 0,ni ! eqn.8 Scherliess draft
a_klnm(i) = S_aM*(a_hf(i)-a_lf(i))/90.+
&2.*a_lf(i)- a_hf(i)
! for testing like in original code
! a_klnm(i)=S_a*(a_hf(i)-a_lf(i))/90.+2.*a_lf(i)-a_hf(i)
! a_klnm(i)=f107d*(a_hf(i)-a_lf(i))/90.+2.*a_lf(i)-a_hf(i)
end do
end subroutine adj_S_a
subroutine constants
!---------------------------------------------------------------
! Purpose: set up constant values (e.g. magnetic grid, convertion
! constants etc)
!
! Method:
!
! Author: A. Maute Nov 2003 am 11/19/03
!--------------------------------------------------------------------
!-------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------
integer :: i,j
real :: fac,lat
rtd = 180./pi ! radians -> deg
dtr = pi/180. ! deg -> radians
sqr2 = sqrt(2.e0)
hr2rd = pi/12. ! pi/12 hrs
dy2rd = 2.*pi/dy2yr ! 2*pi/365.24 average year
deg2mlt = 24./360. ! convert degrees to MLT hours
mlt2deg = 360./24. ! for mlt to mlon
!-------------------------------------------------------------------
! Set grid deltas:
!-------------------------------------------------------------------
dlatm = 180./nmlat
dlonm = 360./nmlon
!-------------------------------------------------------------------
! Set magnetic latitude array
!-------------------------------------------------------------------
do j = 0,nmlat
ylatm(j) = j*dlatm
lat = (ylatm(j) - 90.)*dtr
fac = cos(lat) ! sinIm = 2*sin(lam_m)/sqrt[4-3*cos^2(lam_m)]
fac = 4. - 3.*fac*fac
fac = 2./sqrt( fac )
sinIm_mag(j) = fac*sin( lat )
end do
!------------------------------------------------------------------
! Set magnetic longitude array
!------------------------------------------------------------------
do i = 0,nmlon
ylonm(i) = i*dlonm
end do ! i=1,nmlonp1
!-----------------------------------------------------------------
! find boundary index for weimer
!------------------------------------------------------------------
do j = 0,nmlat
nmlat_wei = j
if( bnd_wei <= ylatm(j) ) then
exit
end if
end do
!-------------------------------------------------------------------
! find latitudinal shift
!-------------------------------------------------------------------
do j = 0,nmlat
ilat_sft = j
if( lat_sft <= ylatm(j) ) then
exit
end if
end do
!------------------------------------------------------------------
! find index for linear interpolation of ed2 at mag.equator
! use 12 deg - same as in TIEGCM
!------------------------------------------------------------------
do j = 0,nmlat
lat = ylatm(j) - 90.
if( lat <= -12. ) then
jmin = j
else if( lat > 12. ) then
jmax = j
exit
end if
end do
end subroutine constants
subroutine prep_fk
!-------------------------------------------------------------------
! Purpose: set up constants factors for f_-k(day) used for empirical model
! to calculate the electric potential
!
! Method:
!
! Author: A. Maute Nov 2003 am 11/19/03
!-------------------------------------------------------------------
ft(1,0) = .75*sqrt( 6.e0 )/pi
ft(1,1) = 2.e0*ft(1,0)
ft(1,2) = 1.e0
ft(2,0) = ft(1,0)
ft(2,1) = -ft(1,1)
ft(2,2) = 1.e0
ft(3,0) = ft(2,1)
ft(3,1) = 0.
ft(3,2) = 1.e0
end subroutine prep_fk
subroutine set_fkflfs( fk, fl, fs )
!------------------------------------------------------------------
! Purpose: set f_-k(day) depending on seasonal flag used for empirical model
! to calculate the electric potential
!
! Method:
!
! Author: A. Maute Nov 2003 am 11/20/03
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------
real, intent(out) ::
& fk(0:2), ! f_-k(day)
& fl(-2:2), ! f_l(ut)
& fs(2) ! f_s(f10.7)
!------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
integer :: lp
real :: ang
real :: lon_ut
!------------------------------------------------------------------
! f_-k(day)
! use factors for iseasav == 0 - Scherliess had iseasav as an input parameter
!------------------------------------------------------------------
lp = iseasav
if( iseasav == 0 ) then
ang = (day + 9.)*dy2rd
fk(0) = sqr2*cos( 2.*ang )
fk(1) = sqr2*cos( ang )
fk(2) = 1.
else if( iseasav >= 1 .and. iseasav <= 3 ) then
fk(0) = ft(lp,0)
fk(1) = ft(lp,1)
fk(2) = ft(lp,2)
else if( iseasav == 4 ) then
fk(0) =0.
fk(1) =0.
fk(2) =1.
end if
!-----------------------------------------------------------------
! f_l(ut)
!-----------------------------------------------------------------
lon_ut = 15.*ut ! 15.*mlt - xmlon + 69.
call ff( lon_ut, 2, fl )
if( iutav ) then ! UT-averaging
ang = fl(0)
fl(:) = 0.
fl(0) = ang
end if
!-----------------------------------------------------------------
! f_s(f10.7) only fs(1) used
!-----------------------------------------------------------------
fs(1) = 1.
! fs(2) = S_a
fs(2) = f107d
end subroutine set_fkflfs
subroutine efield_mid( mlat, mlon, pot )
!------------------------------------------------------------------
! Purpose: calculate the electric potential for low and
! midlatitudes from an empirical model (Scherliess 1999)
!
! Method:
!
! Author: A. Maute Nov 2003 am 11/20/03
!-------------------------------------------------------------------
!------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------
real, intent(in) :: mlat, mlon
real, intent(out) :: pot ! electric potential (V)
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
integer :: i, mp, np, nn
real :: mod_mlat, ct, x
real :: fk(0:2) ! f_-k(day)
real :: fl(-2:2) ! f_l(ut)
real :: fs(2) ! f_s(f10.7)
real :: f(-18:18)
real :: p(0:nm,0:mm) ! P_n^m spherical harmonics
pot = 0. ! initialize
mod_mlat = mlat
if( abs(mlat) <= 0.5 ) then
mod_mlat = 0.5 ! avoid geomag.equator
end if
!------------------------------------------------------------------
! set f_-k, f_l, f_s depending on seasonal flag
!------------------------------------------------------------------
call set_fkflfs( fk, fl, fs )
!------------------------------------------------------------------
! spherical harmonics
!------------------------------------------------------------------
ct = cos( (90. - mod_mlat)*dtr ) ! magnetic colatitude
call pnm( ct, p ) ! calculate P_n^m
call ff( mlon, 18, f ) ! calculate f_m (phi) why 18 if N=12
do i = 0,imax
mp = mf(i)
np = nf(i)
nn = abs(mp) ! P_n^m = P_n^-m
x = a_klnm(i)* fl(lf(i)) * fk(kf(i)) * fs(jf(i))
pot = pot + x*f(mp)*p(np,nn)
end do
end subroutine efield_mid
subroutine prep_weimer
!-----------------------------------------------------------------
! Purpose: for Weimer model calculate IMF angle, IMF magnitude
! tilt of earth
!
! Method: using functions and subroutines from Weimer Model 1996
! output: angle, & ! IMF angle
! bt, & ! IMF magnitude
! tilt ! tilt of earth
!
! Author: A. Maute Nov 2003 am 11/20/03
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
real ::
& angle, ! IMF angle
& bt, ! IMF magnitude
& tilt ! tilt of earth
!-----------------------------------------------------------------
! function declarations
!-----------------------------------------------------------------
real, external :: get_tilt ! in wei96.f
if( by == 0. .and. bz == 0.) then
angle = 0.
else
angle = atan2( by,bz )
end if
angle = angle*rtd
call adjust( angle )
bt = sqrt( by*by + bz*bz )
!-------------------------------------------------------------------
! use month and day of month - calculated with average no.of days per month
! as in Weimer
!-------------------------------------------------------------------
! if(debug) write(iulog,*) 'prep_weimer: day->day of month',
! &iday,imo,iday_m,ut
tilt = get_tilt( iyear, imo, iday_m, ut )
! if(debug) then
! write(iulog,"(/,'efield prep_weimer:')")
! write(iulog,*) ' Bz =',bz
! write(iulog,*) ' By =',by
! write(iulog,*) ' Bt =',bt
! write(iulog,*) ' angle=',angle
! write(iulog,*) ' VSW =',v_sw
! write(iulog,*) ' tilt =',tilt
! end if
call SetModel( angle, bt, tilt, v_sw )
end subroutine prep_weimer
subroutine pot_latsmo( pot, idlat ) ! pots == pot_highlats
!--------------------------------------------------------------------
! Purpose: smoothing in latitude of potential
!
! Method: weighted smoothing in latitude
! assume regular grid spacing
!
! Author: A. Maute Nov 2003 am 11/20/03
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------
integer, intent(in) :: idlat
real, intent(inout) :: pot(0:nmlon,0:nmlat)
!-------------------------------------------------------------------
! local variables
!------------------------------------------------------------------
integer :: ilon, ilat, id
real :: wgt, del
real :: w(-idlat:idlat)
! real :: pot_smo(0:nmlat) ! temp array for smooth. potential
real :: pot_smo(0:nmlon,0:nmlat_wei) ! temp array for smooth. potential
!------------------------------------------------------------------
! weighting factors (regular grid spacing)
!------------------------------------------------------------------
wgt = 0.
do id = -idlat,idlat
del = abs(id)*dlatm ! delta lat_m
w(id) = 1./(del + 1.)
wgt = wgt + w(id)
end do
wgt = 1./wgt
! do ilon = 0,nmlon
! do ilat = idlat,nmlat_wei-idlat
! do ilat = idlat,nmlat-idlat
! pot_smo(ilat) = 0.
! do id = -idlat,idlat ! org. was degree now grid points
! pot_smo(ilat) = pot_smo(ilat) + w(id)*pot(ilon,ilat+id)
! write(iulog,"('pot_latsmo: ilon=',i3,' ilat=',i3,' id=',i3,' pot(ilon,ilat+id)=',e12.4)") ilon,ilat,id,pot(ilon,ilat+id)
! end do
! pot_smo(ilat) = pot_smo(ilat)*wgt
! pot_smo(nmlat-ilat) = pot_smo(ilat)
! end do
! pot(ilon,idlat:nmlat-idlat) = & ! copy back into pot
! pot_smo(idlat:nmlat-idlat)
! pot(ilon,idlat:nmlat_wei-idlat) = pot_smo(idlat:nmlat_wei-idlat)
! pot(ilon,nmlat-nmlat_wei+idlat:nmlat) = pot_smo(nmlat-nmlat_wei+idlat:nmlat)
! pot(ilon,nmlat-nmlat_wei+idlat:nmlat-idlat) = pot_smo(nmlat-nmlat_wei+idlat:nmlat-idlat)
! end do
!$omp parallel do private(ilat)
do ilat = idlat,nmlat_wei-idlat
pot_smo(:,ilat) = matmul( pot(:,ilat-idlat:ilat+idlat),w )*wgt
end do
do ilat = idlat,nmlat_wei-idlat
pot(:,ilat) = pot_smo(:,ilat)
pot(:,nmlat-ilat) = pot_smo(:,ilat)
end do
end subroutine pot_latsmo
subroutine pot_latsmo2( pot, idlat )
!------------------------------------------------------------------
! Purpose: smoothing in latitude of potential
!
! Method: weighted smoothing in latitude
! assume regular grid spacing
!
! Author: A. Maute Nov 2003 am 11/20/03
!------------------------------------------------------------------
!------------------------------------------------------------------
! ... dummy arguments
!------------------------------------------------------------------
integer, intent(in) :: idlat
real, intent(inout) :: pot(0:nmlon,0:nmlat)
!------------------------------------------------------------------
! local variables
!------------------------------------------------------------------
integer :: ilon, ilat, id
real :: wgt, del
real :: w(-idlat:idlat)
! real :: pot_smo(0:nmlat) ! temp array for smooth. potential
real :: pot_smo(0:nmlon,0:nmlath) ! temp array for smooth. potential
!-------------------------------------------------------------------
! weighting factors (regular grid spacing)
!-------------------------------------------------------------------
wgt = 0.
do id = -idlat,idlat
del = abs(id)*dlatm ! delta lat_m
w(id) = 1./(del + 1.)
wgt = wgt + w(id)
end do
wgt = 1./wgt
! do ilon = 0,nmlon
! do ilat = idlat,nmlath-idlat ! ilat = 5:175
! pot_smo(ilat) = 0.
! do id = -idlat,idlat ! org. was degree now grid points
! pot_smo(ilat) = pot_smo(ilat) + w(id)*pot(ilon,ilat+id)
! end do
! pot_smo(ilat) = pot_smo(ilat)*wgt
! end do
! pot(ilon,idlat:nmlath-idlat) = pot_smo(idlat:nmlath-idlat) ! copy back into pot
! end do
!$omp parallel do private(ilat)
do ilat = idlat,nmlath-idlat
pot_smo(:,ilat) = matmul( pot(:,ilat-idlat:ilat+idlat),w )*wgt
end do
do ilat = idlat,nmlath-idlat
pot(:,ilat) = pot_smo(:,ilat)
end do
end subroutine pot_latsmo2
subroutine pot_lonsmo( pot, idlon )
!-------------------------------------------------------------------
! Purpose: smoothing in longitude of potential
!
! Method: weighted smoothing in longitude
! assume regular grid spacing
!
! Author: A. Maute Nov 2003 am 11/20/03
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------
integer, intent(in) :: idlon
real, intent(inout) :: pot(0:nmlon,0:nmlat)
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
integer :: ilon, ilat, id, iabs
real :: wgt, del
real :: w(-idlon:idlon)
real :: pot_smo(0:nmlath) ! temp array for smooth. potential
real :: tmp(-idlon:nmlon+idlon) ! temp array for smooth. potential
!-------------------------------------------------------------------
! weighting factors (regular grid spacing)
!-------------------------------------------------------------------
wgt = 0.
do id = -idlon,idlon
del = abs(id)*dlonm ! delta lon_m
w(id) = 1./(del + 1.)
wgt = wgt + w(id)
end do
wgt = 1./wgt
!-------------------------------------------------------------------
! averaging
!-------------------------------------------------------------------
! do ilon = 0,nmlon
! do ilat = 0,nmlath
! pot_smo(ilat) = 0.
! do id = -idlon,idlon ! org. was degree now grid points
! iabs = ilon + id
! if( iabs > nmlon ) then
! iabs = iabs - nmlon ! test if wrap around
! end if
! if( iabs < 0 ) then
! iabs = iabs + nmlon ! test if wrap around
! end if
! pot_smo(ilat) = pot_smo(ilat) + w(id)*pot(iabs,ilat)
! end do
! pot_smo(ilat) = pot_smo(ilat)*wgt
! pot(ilon,ilat) = pot_smo(ilat) ! copy back into pot
! pot(ilon,nmlat-ilat) = pot_smo(ilat) ! copy back into pot
! end do
! end do
! print*,'www7','pot_efield',pot(149,66),idlon,w
!$omp parallel do private(ilat,ilon,tmp)
do ilat = 0,nmlath
tmp(0:nmlon) = pot(0:nmlon,ilat)
tmp(-idlon:-1) = pot(nmlon-idlon:nmlon-1,ilat)
tmp(nmlon+1:nmlon+idlon) = pot(1:idlon,ilat)
do ilon = 0,nmlon
pot(ilon,ilat)=dot_product(tmp(ilon-idlon:ilon+idlon),w)*wgt
pot(ilon,nmlat-ilat) = pot(ilon,ilat)
! if(ilon.eq.149.and.nmlat-ilat.eq.66)
! &print*,'www9',pot(ilon,nmlat-ilat),tmp(ilon-idlon:ilon+idlon),wgt
end do
end do
! print*,'www8','pot_efield',pot(149,66)
end subroutine pot_lonsmo
subroutine highlat_getbnd( ihlat_bnd )
!------------------------------------------------------------------
! Purpose: calculate the height latitude bounday index ihl_bnd
!
! Method:
! 1. calculate E field from weimar model
! boundary is set where the total electric field exceeds
! 0.015 V/m (corresp. approx. to 300 m/s)
! 2. moved halfways to 54 deg not necessarily equatorwards as in the
! original comment from L. Scherliess- or?
!
! Author: A. Maute Nov 2003 am 11/20/03
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------
integer, intent(out) :: ihlat_bnd(0:nmlon)
!------------------------------------------------------------------
! local variables
!------------------------------------------------------------------
integer :: ilon, ilat, ilat_sft_rvs
real :: mlat, mlt, es, ez, e_tot
ilat_sft_rvs = nmlath - ilat_sft ! pole =0, equ=90
!$omp parallel do private(ilat,ilon,mlt,mlat,es,ez,e_tot)
do ilon = 0,nmlon ! long.
ihlat_bnd(ilon) = 0
mlt = ylonm(ilon)*deg2mlt ! mag.local time ?
do ilat = nmlat_wei+1,0,-1 ! lat. loop moving torwards pole
mlat = 90. - ylatm(ilat) ! mag. latitude pole = 90 equator = 0
call gecmp( mlat, mlt, es, ez ) ! get electric field
e_tot = sqrt( es**2 + ez**2 )
if( abs(e_tot) >= ef_max ) then ! e-filed > limit -> boundary
ihlat_bnd(ilon) = ilat - (ilat - ilat_sft_rvs)/2 ! shift boundary to lat_sft (54deg)
exit
end if
end do
end do
! write(iulog,"('highlat_getbnd: ihlat_bnd=',/,(12i6))") ihlat_bnd
end subroutine highlat_getbnd
subroutine bnd_sinus( ihlat_bnd, itrans_width )
!------------------------------------------------------------------
! Purpose:
! 1. adjust high latitude boundary (ihlat_bnd) sinusoidally
! 2. width of transition zone from midlatitude potential to high latitude
! potential (itrans_width)
!
! Method:
! 1.adjust boundary sinusoidally
! max. wave number to be represented nmax_sin
! RHS(mi) = Sum_phi Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi)*hlat_bnd(phi)
! U(mi,mk) = Sum_phi Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi) *
! Sum_(mk=-nmax_sin)^_(mk=nmax_sin) f_mk(phi)
! single values decomposition of U
! solving U*LSG = RHS
! calculating hlat_bnd:
! hlat_bnd = Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi)*LSG(mi)
!
! 2. width of transition zone from midlatitude potential to high latitude
! potential
! trans_width(phi)=8.-2.*cos(phi)
!
! Author: A. Maute Nov 2003 am 11/20/03
!------------------------------------------------------------------
! use sv_decomp, only : svdcmp, svbksb
!----------------------------------------------------------------------------
! ... dummy arguments
!----------------------------------------------------------------------------
integer, intent(inout) :: ihlat_bnd(0:nmlon) ! loaction of boundary
integer, intent(out) :: itrans_width(0:nmlon) ! width of transition zone
!-----------------------------------------------------------------
! local variables
!-----------------------------------------------------------------
integer, parameter :: nmax_a = 2*nmax_sin+1 ! absolute array length
integer, parameter :: ishf = nmax_sin+1
integer :: ilon, i, i1, j, bnd
real :: sum, mlon
real :: rhs(nmax_a)
real :: lsg(nmax_a)
real :: u(nmax_a,nmax_a)
real :: v(nmax_a,nmax_a)
real :: w(nmax_a,nmax_a)
real :: f(-nmax_sin:nmax_sin,0:nmlon)
!------------------------------------------------------------------
! Sinusoidal Boundary calculation
!------------------------------------------------------------------
rhs(:) = 0.
lsg(:) = 0.
u(:,:) = 0.
v(:,:) = 0.
w(:,:) = 0.
do ilon = 0,nmlon ! long.
bnd = nmlath - ihlat_bnd(ilon) ! switch from pole=0 to pole =90
call ff( ylonm(ilon), nmax_sin, f(-nmax_sin,ilon) )
do i = -nmax_sin,nmax_sin
i1 = i + ishf
rhs(i1) = rhs(i1) + f(i,ilon) * bnd
! write(iulog,*) 'rhs ',ilon,i1,bnd,f(i,ilon),rhs(i+ishf)
do j = -nmax_sin,nmax_sin
u(i1,j+ishf) = u(i1,j+ishf) + f(i,ilon)*f(j,ilon)
! write(iulog,*) 'u ',ilon,i1,j+ishf,u(i+ishf,j+ishf)
end do
end do
end do
! if (debug) write(iulog,*) ' Single Value Decomposition'
call svdcmp( u, nmax_a, nmax_a, nmax_a, nmax_a, w, v )
! if (debug) write(iulog,*) ' Solving'
call svbksb( u, w, v, nmax_a, nmax_a, nmax_a, nmax_a, rhs, lsg )
!
do ilon = 0,nmlon ! long.
! sum = 0.
sum = dot_product( lsg(-nmax_sin+ishf:nmax_sin+ishf),
&f(-nmax_sin:nmax_sin,ilon) )
! do i = -nmax_sin,nmax_sin
! sum = sum + lsg(i+ishf)*f(i,ilon)
! end do
ihlat_bnd(ilon) = nmlath - int( sum + .5 ) ! closest point
itrans_width(ilon) =
&int( 8. - 2.*cos( ylonm(ilon)*dtr ) + .5 )/dlatm ! 6 to 10 deg.
end do
! write(iulog,"('bnd_sinus: ihlat_bnd=',/,(12i6))") ihlat_bnd
! write(iulog,"('bnd_sinus: itrans_width=',/,(12i6))") itrans_width
end subroutine bnd_sinus
subroutine highlat_adjust( pot_highlats, pot_highlat,
&pot_midlat, ihlat_bnd )
!------------------------------------------------------------------
! Purpose: Adjust mid/low latitude electric potential and high latitude
! potential such that there are continous across the mid to high
! latitude boundary
!
! Method:
! 1. integrate Phi_low/mid(phi,bnd) along the boundary mid to high latitude
! 2. integrate Phi_high(phi,bnd) along the boundary mid to high latitude
! 3. adjust Phi_high by delta =
! Int_phi Phi_high(phi,bnd) d phi/360. - Int_phi Phi_low/mid(phi,bnd) d phi/360.
!
! Author: A. Maute Nov 2003 am 11/21/03
!------------------------------------------------------------------
!------------------------------------------------------------------
! ... dummy arguments
!------------------------------------------------------------------
integer, intent(in) :: ihlat_bnd(0:nmlon) ! boundary mid to high latitude
real, intent(in) :: pot_midlat(0:nmlon,0:nmlat) ! low/mid latitude potentail
real, intent(inout) :: pot_highlat(0:nmlon,0:nmlat)! high_lat potential
real, intent(inout) :: pot_highlats(0:nmlon,0:nmlat)! high_lat potential! smoothed high_lat potential
!------------------------------------------------------------------
! local:
!------------------------------------------------------------------
integer :: bnd, ilon, ilat, ilatS, ibnd60, ibnd_hl
real :: pot60, pot_hl, del
!-------------------------------------------------------------------
! 1. integrate Phi_low/mid(phi,bnd) along the boundary mid to high latitude
! 2. integrate Phi_high(phi,bnd) along the boundary mid to high latitude
!-------------------------------------------------------------------
pot60 = 0.
pot_hl = 0.
do ilon = 1,nmlon ! long. ! bnd -> eq to pole 0:90
ibnd60 = nmlat - ihlat_bnd(ilon) ! 0:180 pole to pole
ibnd_hl = ihlat_bnd(ilon) ! colatitude
pot60 = pot60 + pot_midlat(ilon,ibnd60)
pot_hl = pot_hl + pot_highlats(ilon,ibnd_hl)
end do
pot60 = pot60/(nmlon)
pot_hl = pot_hl/(nmlon)
! print*,'www300',pot60,pot_hl,nmlat_wei,nmlon
! if (debug) write(iulog,*) 'Mid-Latitude Boundary Potential =',
! &pot60
! if (debug) write(iulog,*) 'High-Latitude Boundary Potential=',
! &pot_hl
!-------------------------------------------------------------------
! 3. adjust Phi_high by delta =
! Int_phi Phi_high(phi,bnd) d phi/360. - Int_phi Phi_low/mid(phi,bnd) d phi/360.
!-------------------------------------------------------------------
del = pot_hl - pot60
!$omp parallel do private(ilat,ilon,ilats)
do ilat = 0,nmlat_wei ! colatitude
ilats = nmlat - ilat
do ilon = 0,nmlon
pot_highlat(ilon,ilat) = pot_highlat(ilon,ilat) - del
pot_highlat(ilon,ilats) = pot_highlat(ilon,ilats) - del
pot_highlats(ilon,ilat) = pot_highlats(ilon,ilat) - del
pot_highlats(ilon,ilats) = pot_highlats(ilon,ilats) - del
end do
end do
end subroutine highlat_adjust
subroutine interp_poten( pot_highlats, pot_highlat, pot_midlat,
&ihlat_bnd, itrans_width )
!-------------------------------------------------------------------
! Purpose: construct a smooth global electric potential field
!
! Method: construct one global potential field
! 1. low/mid latitude: |lam| < bnd-trans_width
! Phi(phi,lam) = Phi_low(phi,lam)
! 2. high latitude: |lam| > bnd+trans_width
! Phi(phi,lam) = Phi_hl(phi,lam)
! 3. transition zone: bnd-trans_width <= lam <= bnd+trans_width
! a. interpolate between high and low/midlatitude potential
! Phi*(phi,lam) = 1/15*[ 5/(2*trans_width) * {Phi_low(phi,bnd-trans_width)*
! [-lam+bnd+trans_width] + Phi_hl(phi,bnd+trans_width)*
! [lam-bnd+trans_width]} + 10/(2*trans_width) {Phi_low(phi,lam)*
! [-lam+bnd+trans_width] + Phi_hl(phi,lam)*
! [lam-bnd+trans_width]}]
! b. Interpolate between just calculated Potential and the high latitude
! potential in a 3 degree zone poleward of the boundary:
! bnd+trans_width < lam <= bnd+trans_width+ 3 deg
! Phi(phi,lam) = 1/3 { [3-(lam-bnd-trans_width)]* Phi*(phi,lam) +
! [lam-bnd-trans_width)]* Phi_hl*(phi,lam) }
!
! Author: A. Maute Nov 2003 am 11/21/03
!------------------------------------------------------------------
!------------------------------------------------------------------
! ... dummy arguments
!------------------------------------------------------------------
integer, intent(in) :: ihlat_bnd(0:nmlon)
integer, intent(in) :: itrans_width(0:nmlon)
real, intent(in) :: pot_highlats(0:nmlon,0:nmlat)
real, intent(in) :: pot_highlat(0:nmlon,0:nmlat)
real, intent(in) :: pot_midlat(0:nmlon,0:nmlat)
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
real, parameter :: fac = 1./3.
integer :: ilon, ilat
integer :: ibnd, tw, hb1, hb2, lat_ind
integer :: j1, j2
real :: a, b, lat, b1, b2
real :: wrk1, wrk2
!$omp parallel do private(ilat,ilon,ibnd,tw)
do ilon = 0,nmlon
ibnd = ihlat_bnd(ilon) ! high latitude boundary index
tw = itrans_width(ilon) ! width of transition zone (index)
!-------------------------------------------------------------------
! 1. low/mid latitude: |lam| < bnd-trans_width
! Phi(phi,lam) = Phi_low(phi,lam)
!-------------------------------------------------------------------
do ilat = 0,nmlath-(ibnd+tw+1)
potent(ilon,nmlath+ilat) = pot_midlat(ilon,nmlath+ilat)
potent(ilon,nmlath-ilat) = pot_midlat(ilon,nmlath+ilat)
end do
!------------------------------------------------------------------
! 2. high latitude: |lam| > bnd+trans_width
! Phi(phi,lam) = Phi_hl(phi,lam)
!------------------------------------------------------------------
do ilat = 0,ibnd-tw-1
potent(ilon,ilat) = pot_highlats(ilon,nmlat-ilat)
potent(ilon,nmlat-ilat) = pot_highlats(ilon,nmlat-ilat)
end do
end do
!------------------------------------------------------------------
! 3. transition zone: bnd-trans_width <= lam <= bnd+trans_width
!------------------------------------------------------------------
! a. interpolate between high and low/midlatitude potential
! update only southern hemisphere (northern hemisphere is copied
! after smoothing)
!------------------------------------------------------------------
!!$omp parallel do private(ilat,ilon,ibnd,tw,a,b,b1,b2,hb1,hb2,
! &lat_ind,j1,j2,wrk1,wrk2)
do ilon = 0,nmlon
ibnd = ihlat_bnd(ilon) ! high latitude boundary index
tw = itrans_width(ilon) ! width of transition zone (index)
a = 1./(2.*tw)
b1 = (nmlath - ibnd + tw)*a
b2 = (nmlath - ibnd - tw)*a
hb1 = nmlath - (ibnd + tw)
j1 = nmlath - hb1
hb2 = nmlath - (ibnd - tw)
j2 = nmlath - hb2
wrk1 = pot_midlat(ilon,j1)
wrk2 = pot_highlats(ilon,j2)
! write(iulog,*) 'pot_all ',ilon,hb1,hb2,nmlath -ibnd,tw
do ilat = ibnd-tw,ibnd+tw
lat_ind = nmlath - ilat
potent(ilon,ilat) =
& fac*((wrk1 + 2.*pot_midlat(ilon,ilat))*(b1 - a*lat_ind)
& + (wrk2 + 2.*pot_highlats(ilon,ilat))*(a*lat_ind - b2))
potent(ilon,nmlat-ilat) = potent(ilon,ilat)
end do
!------------------------------------------------------------------
! b. Interpolate between just calculated Potential and the high latitude
! potential in a 3 degree zone poleward of the boundary
!------------------------------------------------------------------
do ilat = hb2+1,nmlath
a = max( 3./dlatm - (ilat - hb2 - 1),0. )
b = 3./dlatm - a
potent(ilon,nmlath-ilat) = (a*potent(ilon,nmlath-ilat)
& + b*pot_highlat(ilon,nmlath-ilat))/3.*dlatm
potent(ilon,nmlath+ilat) = potent(ilon,nmlath-ilat)
end do
end do
end subroutine interp_poten
subroutine DerivPotential
!-----------------------------------------------------------------
! Purpose: calulates the electric field [V/m] from the electric potential
!
! Method: Richmond [1995] eqn 5.9-5.10
! ed1(:,:) = Ed1 = - 1/[R cos lam_m] d PHI/d phi_m
! ed2(:,:) = Ed2 = 1/R d PHI/d lam_m /sin I_m
! R = R_e + h_r we assume a reference height of 130 km which is also
! used in the TIEGCM code
!
! Author: A. Maute Dec 2003 am 12/16/03
!-----------------------------------------------------------------
integer :: i, j, ip1f, ip2f, ip3f
real :: coslm, r, fac, wrk
real :: wrk1d(0:nmlon)
r = r_e + h_r ! earth radius + reference height [m]
!-----------------------------------------------------------------
! ed2= Ed2 is the equatorward/downward component of the electric field, at all
! geomagnetic grid points (central differencing)
!-----------------------------------------------------------------
fac = .5/(dlatm*dtr*r)
!$omp parallel do private(j, i, wrk )
do j = 1,nmlath-1 ! southern hemisphere
! idea
wrk = fac/sinIm_mag(j)
! wrk = fac
do i = 0,nmlon
ed2(i,j) = (potent(i,j+1) - potent(i,j-1))*wrk
end do
end do
!$omp parallel do private(j, i, wrk )
do j = nmlath+1,nmlat-1 ! northern hemisphere
wrk = fac/sinIm_mag(j)
do i = 0,nmlon
ed2(i,j) = (potent(i,j+1) - potent(i,j-1))*wrk
end do
end do
!-----------------------------------------------------------------------
! Interpolate of ed2 between between -12 <= lam_m <= 12 degrees:
!-----------------------------------------------------------------------
wrk1d(:) = ed2(:,jmax) - ed2(:,jmin)
do j = jmin+1,jmax-1
fac = (ylatm(j) - ylatm(jmin))/(ylatm(jmax) - ylatm(jmin))
do i = 0,nmlon
ed2(i,j) = ed2(i,jmin) + fac*wrk1d(i)
end do
end do
!-----------------------------------------------------------------------
! ed1= Ed1 is the zonal component of the electric field, at all
! geomagnetic grid points (central differencing)
!-----------------------------------------------------------------------
fac = .5/(dlonm*dtr*r)
!$omp parallel do private(j, i, wrk, coslm )
do j = 1,nmlat-1
coslm = ylatm(j) - 90.
coslm = cos( coslm*dtr )
wrk = fac/coslm
do i = 1,nmlon-1
ed1(i,j) = -(potent(i+1,j) - potent(i-1,j))*wrk
end do
i = 0
ed1(i,j) = -(potent(i+1,j) - potent(nmlon-1,j))*wrk
ed1(nmlon,j) = ed1(i,j)
end do
!-----------------------------------------------------------------------
! Poles:
!-----------------------------------------------------------------------
do i = 0,nmlon
ip1f = i + nmlon/4
if( ip1f > nmlon ) then
ip1f = ip1f - nmlon
end if
ip2f = i + nmlon/2
if( ip2f > nmlon ) then
ip2f = ip2f - nmlon
end if
ip3f = i + 3*nmlon/4
if( ip3f > nmlon ) then
ip3f = ip3f - nmlon
end if
ed1(i,0)=.25*(ed1(i,1)-ed1(ip2f,1)+ed2(ip1f,1)-ed2(ip3f,1))
ed1(i,nmlat) = .25*(ed1(i,nmlat-1) - ed1(ip2f,nmlat-1)
& + ed2(ip1f,nmlat-1) - ed2(ip3f,nmlat-1))
ed2(i,0)=.25*(ed2(i,1)-ed2(ip2f,1)-ed1(ip1f,1)+ed1(ip3f,1))
ed2(i,nmlat) = .25*(ed2(i,nmlat-1) - ed2(ip2f,nmlat-1)
& - ed1(ip1f,nmlat-1) + ed1(ip3f,nmlat-1))
end do
end subroutine DerivPotential
end module efield
!
! Purpose:
! Subroutines to calculate the electric potentials from the Weimer '96 model of
! the polar cap ionospheric electric potentials.
!
! Method:
!
! To use, first call subroutine ReadCoef once.
! Next, call SetModel with the specified input parameters.
! The function EpotVal(gLAT,gMLT) can then be used repeatively to get the
! electric potential at the desired location in geomagnetic coordinates.
! Subroutines to calculate the electric potentials from the Weimer '96 model of
! the polar cap ionospheric electric potentials.
!
!
! Author: A. Maute Dec 2003
! This code is protected by copyright and is
! distributed for research or educational use only.
! Commerical use without written permission from Dan Weimer/MRC is prohibited.
!
!*********************** Copyright 1996, Dan Weimer/MRC ***********************
!==================================================================
FUNCTION EpotVal(gLAT,gMLT)
!
!-----------------------------------------------------------------------
! Return the value of the electric potential in kV at
! corrected geomagnetic coordinates gLAT (degrees) and gMLT (hours).
!
! Must first call ReadCoef and SetModel to set up the model coeficients for
! the desired values of Bt, IMF clock angle, Dipole tilt angle, and SW Vel.
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only: Coef =>Coef,ML=>ML,MM=>MM1
implicit none
!
!-----------------------------Return Value------------------------------
!
real EpotVal
!
!-------------------------------Commons---------------------------------
!
! INTEGER ML,MM
! REAL Coef(0:1,0:8,0:3),pi
! COMMON/SetCoef/Coef,pi,ML,MM
real pi
!
!------------------------------Arguments--------------------------------
!
REAL gLAT,gMLT
!
!---------------------------Local variables-----------------------------
!
integer limit,l,m
Real Theta,Phi,Z,ct,Phim
real r
REAL Plm(0:20,0:20)
!
!-----------------------------------------------------------------------
!
pi=3.141592653
r=90.-gLAT
IF(r .LT. 45.)THEN
Theta=r*pi/45.
Phi=gMLT*pi/12.
Z=Coef(0,0,0)
ct=COS(Theta)
CALL Legendre(ct,ML,MM,Plm)
DO l=1,ML
Z=Z + Coef(0,l,0)*Plm(l,0)
IF(l.LT.MM)THEN
limit=l
ELSE
limit=MM
ENDIF
DO m=1,limit
phim=phi*m
Z=Z + Coef(0,l,m)*Plm(l,m)*COS(phim) +
& Coef(1,l,m)*Plm(l,m)*SIN(phim)
ENDDO
ENDDO
ELSE
Z=0.
ENDIF
EpotVal=Z
! print*,'www0',Z,Coef,Plm,ct
RETURN
END FUNCTION EpotVal
!================================================================================================
SUBROUTINE ReadCoef (wei96_file)
!
!-----------------------------------------------------------------------
!
! Read in the data file with the model coefficients
!
!*********************** Copyright 1996, Dan Weimer/MRC ***********************
!
! NCAR addition (Jan 97): initialize constants used in GECMP
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
! use ioFileMod, only : getfil
! use units, only : getunit, freeunit
! use abortutils, only : endrun
! use cam_logfile, only : iulog
use efield, only: ALAMN =>ALAMN,ALAMX=>ALAMX,ALAMR=>ALAMR,
&STPD=>STPD,STP2=>STP2,CSTP=>CSTP,SSTP=>SSTP,
&Cn=>Cn,MaxL=>MaxL,MaxM=>MaxM,MaxN=>MaxN
implicit none
!
!-------------------------------Commons---------------------------------
!
! real alamn, alamx, alamr, stpd, stp2, cstp, sstp
! COMMON /CECMP/ ALAMN,ALAMX,ALAMR,STPD,STP2,CSTP,SSTP
! ALAMN = Absolute min latitude (deg) of model
! ALAMX = Absolute max latitude (deg) for normal gradient calc.
! STPD = Angular dist (deg) of step @ 300km above earth (r=6371km)
! STP2 = Denominator in gradient calc
!
!------------------------------Arguments--------------------------------
!
character(len=*), intent(in) :: wei96_file
!
!-----------------------------Parameters------------------------------
!
real d2r, r2d
PARAMETER ( D2R = 0.0174532925199432957692369076847 ,
& R2D = 57.2957795130823208767981548147)
!
!---------------------------Local variables-----------------------------
!
INTEGER udat,unit,ios
integer ll,mm,k,m,klimit,kk,nn,ii,i,n,ilimit,mlimit,l
REAL C(0:3)
real stpr, step
CHARACTER*15 skip
INTEGER iulog
! INTEGER MaxL,MaxM,MaxN,iulog
! REAL Cn( 0:3 , 0:1 , 0:4 , 0:1 , 0:8 , 0:3 )
! COMMON /AllCoefs/Cn,MaxL,MaxM,MaxN
character(len=256) :: locfn
!
!-----------------------------------------------------------------------
iulog=14
STEP = 10.
STPR = STEP/6671.
STPD = STPR*R2D
STP2 = 2.*STEP
CSTP = COS (STPR)
SSTP = SQRT (1. - CSTP*CSTP)
ALAMN = 45.
ALAMX = 90. - STPD
ALAMR = ALAMN*D2R
! End NCAR addition
!
! get coeff_file
! unit= getunit()
unit= 600
! print*, 'Weimer: getting file ',trim(wei96_file),
! &' unit ',unit
! call getfil( wei96_file, locfn, 0 )
locfn= wei96_file
!
! write(iulog,*) 'Weimer: opening file ',trim(locfn),
! &' unit ',unit
! OPEN(unit=unit,file=trim(locfn),
open(unit=unit,file=locfn,status = 'old',iostat = ios)
if(ios.gt.0) then
print*, 'Weimer: error in opening wei96.cofcnts',
&' unit ',unit
! call endrun
endif
900 FORMAT(A15)
c1000 FORMAT(3I8)
1000 format(3i8)
2000 FORMAT(3I2)
3000 FORMAT(2I2,4E15.6)
! READ(udat,900) skip
! write(iulog,*) 'Weimer: reading file ',trim(locfn),
! &' unit ',unit
! READ(unit,1000,iostat = ios) MaxL,MaxM,MaxN
read(unit,1000,iostat = ios) MaxL,MaxM,MaxN
! print*,'www0',ios,MaxL,MaxM,MaxN
! if(ios.gt.0) then
! write(iulog,*)
! &'ReadCoef: error in reading wei96.cofcnts file',
! &' unit ',unit
! call endrun
! endif
DO l=0,MaxL
IF(l.LT.MaxM)THEN
mlimit=l
ELSE
mlimit=MaxM
ENDIF
DO m=0,mlimit
IF(m.LT.1)THEN
klimit=0
ELSE
klimit=1
ENDIF
DO k=0,klimit
read(unit,2000,iostat = ios) ll,mm,kk
! print*,k,ll,mm,kk
! if(ios.gt.0) then
! write(iulog,*)
! &'ReadCoef: error in reading wei96.cofcnts file',' unit ',
! &unit
! call endrun
! endif
! IF(ll.NE.l .OR. mm.NE.m .OR. kk.NE.k)THEN
! WRITE(IULOG,*)'Data File Format Error'
! CALL ENDRUN
! ENDIF
DO n=0,MaxN
IF(n.LT.1)THEN
ilimit=0
ELSE
ilimit=1
ENDIF
DO i=0,ilimit
READ(unit,3000,iostat = ios) nn,ii,C
! print*,'www0',nn,ii,C,i,n,k,l,m
! if(ios.gt.0) then
! write(iulog,*) 'ReadCoef: error in reading',
! & ' wei96.cofcnts file',' unit ',unit
! call endrun
! endif
! IF(nn.NE.n .OR. ii.NE.i)THEN
! WRITE(IULOG,*)'Data File Format Error'
! CALL ENDRUN
! ENDIF
Cn(0,i,n,k,l,m)=C(0)
Cn(1,i,n,k,l,m)=C(1)
Cn(2,i,n,k,l,m)=C(2)
Cn(3,i,n,k,l,m)=C(3)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
close(unit)
! call freeunit(unit)
!
RETURN
END SUBROUTINE ReadCoef
!================================================================================================
FUNCTION FSVal(omega,MaxN,FSC)
!
!-----------------------------------------------------------------------
! Evaluate a Sine/Cosine Fourier series for N terms up to MaxN
! at angle omega, given the coefficients in FSC
!
!*********************** Copyright 1996, Dan Weimer/MRC ***************
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!-----------------------------Return Value------------------------------
!
real FSVal
!
!------------------------------Arguments--------------------------------
!
INTEGER MaxN
! REAL omega,FSC(0:1,0:*)
REAL omega,FSC(0:1,0:4)
!
!---------------------------Local variables-----------------------------
!
INTEGER n
REAL Y,theta
!
!-----------------------------------------------------------------------
!
Y=0.
DO n=0,MaxN
theta=omega*n
Y=Y + FSC(0,n)*COS(theta) + FSC(1,n)*SIN(theta)
ENDDO
FSVal=Y
! print*,'www00',Y,FSC
RETURN
END FUNCTION FSVal
!================================================================================================
SUBROUTINE SetModel(angle,Bt,Tilt,SWVel)
!
!-----------------------------------------------------------------------
! Calculate the complete set of spherical harmonic coefficients,
! given an arbitrary IMF angle (degrees from northward toward +Y),
! magnitude Bt (nT), dipole tilt angle (degrees),
! and solar wind velocity (km/sec).
! Returns the Coef in the common block SetCoef.
!
!*********************** Copyright 1996, Dan Weimer/MRC ***********************
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only: Cn=>Cn,MaxL=>MaxL,MaxM=>MaxM,MaxN=>MaxN
&,Coef=>coef,ML=>ML,MM=>MM1
implicit none
!
!-------------------------------Commons---------------------------------
!
! INTEGER MaxL,MaxM,MaxN
! REAL Cn( 0:3 , 0:1 , 0:4 , 0:1 , 0:8 , 0:3 )
! COMMON /AllCoefs/Cn,MaxL,MaxM,MaxN
! INTEGER ML,MM
! REAL Coef(0:1,0:8,0:3),pi
! COMMON/SetCoef/Coef,pi,ML,MM
real pi
!
!------------------------------Arguments--------------------------------
!
REAL angle,Bt,Tilt,SWVel
!
!---------------------------Local variables-----------------------------
!
integer n, k, ilimit, i, klimit, l, m, mlimit
REAL FSC(0:1,0:4), fsval, omega, sintilt
!
!-----------------------------------------------------------------------
!
pi=3.141592653
ML=MaxL
MM=MaxM
SinTilt=SIN(Tilt*pi/180.)
! SinTilt=SIND(Tilt)
omega=angle*pi/180.
fsc(1,0) = 0.
DO l=0,MaxL
IF(l.LT.MaxM)THEN
mlimit=l
ELSE
mlimit=MaxM
ENDIF
DO m=0,mlimit
IF(m.LT.1)THEN
klimit=0
ELSE
klimit=1
ENDIF
DO k=0,klimit
! Retrieve the regression coefficients and evaluate the function
! as a function of Bt,Tilt,and SWVel to get each Fourier coefficient.
DO n=0,MaxN
IF(n.LT.1)THEN
ilimit=0
ELSE
ilimit=1
ENDIF
DO i=0,ilimit
FSC(i,n)=Cn(0,i,n,k,l,m) + Bt*Cn(1,i,n,k,l,m) +
& SinTilt*Cn(2,i,n,k,l,m) + SWVel*Cn(3,i,n,k,l,m)
ENDDO
ENDDO
! Next evaluate the Fourier series as a function of angle.
Coef(k,l,m)=FSVal(omega,MaxN,FSC)
ENDDO
ENDDO
ENDDO
! print*,'www000',FSC(0,0),Cn,Bt,SinTilt,SWVel
RETURN
END SUBROUTINE SetModel
!================================================================================================
SUBROUTINE LEGENDRE(x,lmax,mmax,Plm)
!
!-----------------------------------------------------------------------
! compute Associate Legendre Function P_l^m(x)
! for all l up to lmax and all m up to mmax.
! returns results in array Plm
! if X is out of range ( abs(x)>1 ) then value is returned as if x=1.
!
!*********************** Copyright 1996, Dan Weimer/MRC ***********************
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
! use cam_logfile, only : iulog
implicit none
!
!------------------------------Arguments--------------------------------
!
integer lmax, mmax
real x, Plm(0:20,0:20)
!
!---------------------------Local variables-----------------------------
!
integer m, lm2, l, iulog
real xx, fact
iulog=14
!
!-----------------------------------------------------------------------
!
DO l=0,20
DO m=0,20
Plm(l,m)=0.
ENDDO
ENDDO
xx=MIN(x,1.)
xx=MAX(xx,-1.)
! IF(lmax .LT. 0 .OR. mmax .LT. 0 .OR. mmax .GT. lmax )THEN
! write(iulog,*)'Bad arguments to Legendre'
! RETURN
! ENDIF
! First calculate all Pl0 for l=0 to l
Plm(0,0)=1.
IF(lmax.GT.0)Plm(1,0)=xx
IF (lmax .GT. 1 )THEN
DO L=2,lmax
Plm(L,0)=( (2.*L-1)*xx*Plm(L-1,0) -
&(L-1)*Plm(L-2,0) )/L
ENDDO
ENDIF
IF (mmax .EQ. 0 )RETURN
fact=SQRT( (1.-xx)*(1.+xx) )
DO M=1,mmax
DO L=m,lmax
lm2=MAX(L-2,0)
Plm(L,M)=Plm(lm2,M) - ( 2*L-1)*fact*Plm(L-1,M-1)
ENDDO
ENDDO
RETURN
END SUBROUTINE LEGENDRE
!================================================================================================
!*********************** Copyright 1996, Dan Weimer/MRC ***********************
!CC NCAR MODIFIED (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! The following routines (translib.for) were added to return the dipole tilt. C
! GET_TILT was initially a procedure (TRANS), here it has been changed into C
! a function which returns the dipole tilt. C
! Barbara Emery (emery@ncar.ucar.edu) and William Golesorkhi, HAO/NCAR (3/96) C
!CC NCAR MODIFIED (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! COORDINATE TRANSFORMATION UTILITIES
!**********************************************************************
FUNCTION GET_TILT(YEAR,MONTH,DAY,HOUR)
!
!-----------------------------------------------------------------------
!CC NCAR MODIFIED (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! The following line initially was: C
! SUBROUTINE TRANS(YEAR,MONTH,DAY,HOUR,IDBUG) C
! It has been changed to return the dipole tilt from this function call. C
!CC NCAR MODIFIED (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! THIS SUBROUTINE DERIVES THE ROTATION MATRICES AM(I,J,K) FOR 11
! TRANSFORMATIONS, IDENTIFIED BY K.
! K=1 TRANSFORMS GSE to GEO
! K=2 " GEO to MAG
! K=3 " GSE to MAG
! K=4 " GSE to GSM
! K=5 " GEO to GSM
! K=6 " GSM to MAG
! K=7 " GSE to GEI
! K=8 " GEI to GEO
! K=9 " GSM to SM
! K=10 " GEO to SM
! K=11 " MAG to SM
!
! IF IDBUG IS NOT 0, THEN OUTPUTS DIAGNOSTIC INFORMATION TO
! FILE UNIT=IDBUG
!
! The formal names of the coordinate systems are:
! GSE - Geocentric Solar Ecliptic
! GEO - Geographic
! MAG - Geomagnetic
! GSM - Geocentric Solar Magnetospheric
! SM - Solar Magnetic
!
! THE ARRAY CX(I) ENCODES VARIOUS ANGLES, STORED IN DEGREES
! ST(I) AND CT(I) ARE SINES & COSINES.
!
! Program author: D. R. Weimer
!
! Some of this code has been copied from subroutines which had been
! obtained from D. Stern, NASA/GSFC. Other formulas are from "Space
! Physics Coordinate Transformations: A User Guide" by M. Hapgood (1991).
!
! The formulas for the calculation of Greenwich mean sidereal time (GMST)
! and the sun's location are from "Almanac for Computers 1990",
! U.S. Naval Observatory.
!
!-----------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only:CX=>CX,ST=>ST,CT=>CT,AM=>AM
&,EPOCH=>EPOCH,TH0=>TH0,PH0=>PH0,DIPOLE=>DIPOLE
implicit none
!
!-----------------------------Return Value--------------------------
!
real get_tilt
!
!-------------------------------Commons---------------------------------
!
! real cx, st, ct, am
! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
! real epoch, th0, ph0, dipole
! COMMON/MFIELD/EPOCH,TH0,PH0,DIPOLE
! DATA EPOCH,TH0,PH0,DIPOLE/1980.,11.19,-70.76,.30574/
!
!------------------------------Arguments--------------------------------
!
INTEGER YEAR, MONTH, DAY
REAL HOUR
!
!-----------------------------Parameters------------------------------
!
INTEGER GSEGEO,GEOGSE,GEOMAG,MAGGEO
INTEGER GSEMAG,MAGGSE,GSEGSM,GSMGSE
INTEGER GEOGSM,GSMGEO,GSMMAG,MAGGSM
INTEGER GSEGEI,GEIGSE,GEIGEO,GEOGEI
INTEGER GSMSM,SMGSM,GEOSM,SMGEO,MAGSM,SMMAG
PARAMETER (GSEGEO= 1,GEOGSE=-1,GEOMAG= 2,MAGGEO=-2)
PARAMETER (GSEMAG= 3,MAGGSE=-3,GSEGSM= 4,GSMGSE=-4)
PARAMETER (GEOGSM= 5,GSMGEO=-5,GSMMAG= 6,MAGGSM=-6)
PARAMETER (GSEGEI= 7,GEIGSE=-7,GEIGEO= 8,GEOGEI=-8)
PARAMETER (GSMSM = 9,SMGSM =-9,GEOSM =10,SMGEO=-10)
PARAMETER (MAGSM =11,SMMAG =-11)
!
!---------------------------Local variables-----------------------------
!
integer IDBUG
integer j, k, jd, iyr, i, mjd
REAL UT, T0, GMSTD, GMSTH, ECLIP, MA, LAMD, SUNLON, pi
real b32, b33, b3
!
!-------------------------External Functions----------------------------
!
integer julday_wam
external julday_wam
!
!-----------------------------------------------------------------------
!
! EPOCH=1980.
! TH0=11.19
! PH0=-70.76
! DIPOLE=.30574
pi=3.141592653
! pi=2.*ASIN(1.)
!CC NCAR MODIFICATION (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! IDBUG=0 to prevent printing data to the screen or writing data to a file. C
IDBUG = 0
!CC NCAR MODIFICATION (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IF(YEAR.LT.1900)THEN
IYR=1900+YEAR
ELSE
IYR=YEAR
ENDIF
UT=HOUR
JD=JULDAY_WAM(MONTH,DAY,IYR)
MJD=JD-2400001
! T0=(real(MJD,r8)-51544.5)/36525.0
T0=(float(MJD)-51544.5)/36525.0
GMSTD=100.4606184 +36000.770*T0 +3.87933E-4*T0*T0 +
& 15.0410686*UT
CALL ADJUST(GMSTD)
GMSTH=GMSTD*24./360.
ECLIP=23.439 - 0.013*T0
MA=357.528 + 35999.050*T0 + 0.041066678*UT
CALL ADJUST(MA)
LAMD=280.460 + 36000.772*T0 + 0.041068642*UT
CALL ADJUST(LAMD)
SUNLON=LAMD + (1.915-0.0048*T0)*SIN(MA*pi/180.) + 0.020*
& SIN(2.*MA*pi/180.)
CALL ADJUST(SUNLON)
! IF(IDBUG.NE.0)THEN
! WRITE(IDBUG,*) YEAR,MONTH,DAY,HOUR
! WRITE(IDBUG,*) 'MJD=',MJD
! WRITE(IDBUG,*) 'T0=',T0
! WRITE(IDBUG,*) 'GMSTH=',GMSTH
! WRITE(IDBUG,*) 'ECLIPTIC OBLIQUITY=',ECLIP
! WRITE(IDBUG,*) 'MEAN ANOMALY=',MA
! WRITE(IDBUG,*) 'MEAN LONGITUDE=',LAMD
! WRITE(IDBUG,*) 'TRUE LONGITUDE=',SUNLON
! ENDIF
CX(1)= GMSTD
CX(2) = ECLIP
CX(3) = SUNLON
CX(4) = TH0
CX(5) = PH0
! Derived later:
! CX(6) = Dipole tilt angle
! CX(7) = Angle between sun and magnetic pole
! CX(8) = Subsolar point latitude
! CX(9) = Subsolar point longitude
DO I=1,5
ST(I) = SIN(CX(I)*pi/180.)
CT(I) = COS(CX(I)*pi/180.)
ENDDO
!
AM(1,1,GSEGEI) = CT(3)
AM(1,2,GSEGEI) = -ST(3)
AM(1,3,GSEGEI) = 0.
AM(2,1,GSEGEI) = ST(3)*CT(2)
AM(2,2,GSEGEI) = CT(3)*CT(2)
AM(2,3,GSEGEI) = -ST(2)
AM(3,1,GSEGEI) = ST(3)*ST(2)
AM(3,2,GSEGEI) = CT(3)*ST(2)
AM(3,3,GSEGEI) = CT(2)
!
AM(1,1,GEIGEO) = CT(1)
AM(1,2,GEIGEO) = ST(1)
AM(1,3,GEIGEO) = 0.
AM(2,1,GEIGEO) = -ST(1)
AM(2,2,GEIGEO) = CT(1)
AM(2,3,GEIGEO) = 0.
AM(3,1,GEIGEO) = 0.
AM(3,2,GEIGEO) = 0.
AM(3,3,GEIGEO) = 1.
!
DO I=1,3
DO J=1,3
AM(I,J,GSEGEO) = AM(I,1,GEIGEO)*AM(1,J,GSEGEI) +
&AM(I,2,GEIGEO)*AM(2,J,GSEGEI) +AM(I,3,GEIGEO)*AM(3,J,GSEGEI)
ENDDO
ENDDO
!
AM(1,1,GEOMAG) = CT(4)*CT(5)
AM(1,2,GEOMAG) = CT(4)*ST(5)
AM(1,3,GEOMAG) =-ST(4)
AM(2,1,GEOMAG) =-ST(5)
AM(2,2,GEOMAG) = CT(5)
AM(2,3,GEOMAG) = 0.
AM(3,1,GEOMAG) = ST(4)*CT(5)
AM(3,2,GEOMAG) = ST(4)*ST(5)
AM(3,3,GEOMAG) = CT(4)
!
DO I=1,3
DO J=1,3
AM(I,J,GSEMAG) = AM(I,1,GEOMAG)*AM(1,J,GSEGEO) +
&AM(I,2,GEOMAG)*AM(2,J,GSEGEO) +AM(I,3,GEOMAG)*AM(3,J,GSEGEO)
ENDDO
ENDDO
!
B32 = AM(3,2,GSEMAG)
B33 = AM(3,3,GSEMAG)
B3 = SQRT(B32*B32+B33*B33)
IF (B33.LE.0.) B3 = -B3
!
AM(2,2,GSEGSM) = B33/B3
AM(3,3,GSEGSM) = AM(2,2,GSEGSM)
AM(3,2,GSEGSM) = B32/B3
AM(2,3,GSEGSM) =-AM(3,2,GSEGSM)
AM(1,1,GSEGSM) = 1.
AM(1,2,GSEGSM) = 0.
AM(1,3,GSEGSM) = 0.
AM(2,1,GSEGSM) = 0.
AM(3,1,GSEGSM) = 0.
!
DO I=1,3
DO J=1,3
AM(I,J,GEOGSM) = AM(I,1,GSEGSM)*AM(J,1,GSEGEO) +
&AM(I,2,GSEGSM)*AM(J,2,GSEGEO) +
&AM(I,3,GSEGSM)*AM(J,3,GSEGEO)
ENDDO
ENDDO
!
DO I=1,3
DO J=1,3
AM(I,J,GSMMAG) = AM(I,1,GEOMAG)*AM(J,1,GEOGSM) +
&AM(I,2,GEOMAG)*AM(J,2,GEOGSM) +
&AM(I,3,GEOMAG)*AM(J,3,GEOGSM)
ENDDO
ENDDO
!
ST(6) = AM(3,1,GSEMAG)
CT(6) = SQRT(1.-ST(6)*ST(6))
CX(6) = ASIN(ST(6)*pi/180.)
AM(1,1,GSMSM) = CT(6)
AM(1,2,GSMSM) = 0.
AM(1,3,GSMSM) = -ST(6)
AM(2,1,GSMSM) = 0.
AM(2,2,GSMSM) = 1.
AM(2,3,GSMSM) = 0.
AM(3,1,GSMSM) = ST(6)
AM(3,2,GSMSM) = 0.
AM(3,3,GSMSM) = CT(6)
!
DO I=1,3
DO J=1,3
AM(I,J,GEOSM) = AM(I,1,GSMSM)*AM(1,J,GEOGSM) +
&AM(I,2,GSMSM)*AM(2,J,GEOGSM) +
&AM(I,3,GSMSM)*AM(3,J,GEOGSM)
ENDDO
ENDDO
!
DO I=1,3
DO J=1,3
AM(I,J,MAGSM) = AM(I,1,GSMSM)*AM(J,1,GSMMAG) +
& AM(I,2,GSMSM)*AM(J,2,GSMMAG) +
&AM(I,3,GSMSM)*AM(J,3,GSMMAG)
ENDDO
ENDDO
!
CX(7)=ATAN2( AM(2,1,11) , AM(1,1,11) )
CX(7)=CX(7)*180./pi
CX(8)=ASIN( AM(3,1,1)*pi/180. )
CX(9)=ATAN2( AM(2,1,1) , AM(1,1,1) )
CX(9)=CX(9)*180./pi
IF(IDBUG.NE.0)THEN
! WRITE(IDBUG,*) 'Dipole tilt angle=',CX(6)
! WRITE(IDBUG,*) 'Angle between sun and magnetic pole=',
! &CX(7)
! WRITE(IDBUG,*) 'Subsolar point latitude=',CX(8)
! WRITE(IDBUG,*) 'Subsolar point longitude=',CX(9)
DO K=1,11
! WRITE(IDBUG,1001) K
DO I=1,3
! WRITE(IDBUG,1002) (AM(I,J,K),J=1,3)
ENDDO
ENDDO
1001 FORMAT(' ROTATION MATRIX ',I2)
1002 FORMAT(3F9.5)
ENDIF
!CC NCAR MODIFICATION (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! The next line was added to return the dipole tilt from this function call. C
GET_TILT = CX(6)
!CC NCAR MODIFICATION (3/96) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
RETURN
END FUNCTION GET_TILT
!======================================================================
SUBROUTINE ROTATE (X,Y,Z,I)
!
!-----------------------------------------------------------------------
! THIS SUBROUTINE APPLIES TO THE VECTOR (X,Y,Z) THE ITH ROTATION
! MATRIX AM(N,M,I) GENERATED BY SUBROUTINE TRANS
! IF I IS NEGATIVE, THEN THE INVERSE ROTATION IS APPLIED
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!------------------------------Arguments--------------------------------
!
integer i
REAL X,Y,Z
!
!---------------------------Local variables-----------------------------
!
REAL A(3)
!
!-----------------------------------------------------------------------
!
A(1)=X
A(2)=Y
A(3)=Z
CALL ROTATEV(A,A,I)
X=A(1)
Y=A(2)
Z=A(3)
RETURN
END SUBROUTINE ROTATE
!======================================================================
SUBROUTINE ROTATEV (A,B,I)
!
!-----------------------------------------------------------------------
! THIS SUBROUTINE APPLIES TO THE VECTOR A(3) THE ITH ROTATION
! MATRIX AM(N,M,I) GENERATED BY SUBROUTINE TRANS
! AND OUTPUTS THE CONVERTED VECTOR B(3), WITH NO CHANGE TO A.
! IF I IS NEGATIVE, THEN THE INVERSE ROTATION IS APPLIED
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
! use cam_logfile, only : iulog
! use abortutils, only : endrun
use efield, only:CX=>CX,ST=>ST,CT=>CT,AM=>AM
implicit none
!
!-------------------------------Commons---------------------------------
!
! real cx, st, ct, am
! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
!
!------------------------------Arguments--------------------------------
!
integer i
REAL A(3),B(3)
!
!---------------------------Local variables-----------------------------
!
integer id, j, iulog
real xa, ya, za
iulog=14
!
!-----------------------------------------------------------------------
!
! IF(I.EQ.0 .OR. IABS(I).GT.11)THEN
! WRITE(IULOG,*)'ROTATEV CALLED WITH UNDEFINED TRANSFORMATION'
! CALL ENDRUN
! ENDIF
XA = A(1)
YA = A(2)
ZA = A(3)
IF(I.GT.0)THEN
ID=I
DO J=1,3
B(J) = XA*AM(J,1,ID) + YA*AM(J,2,ID) + ZA*AM(J,3,ID)
ENDDO
ELSE
ID=-I
DO J=1,3
B(J) = XA*AM(1,J,ID) + YA*AM(2,J,ID) + ZA*AM(3,J,ID)
ENDDO
ENDIF
RETURN
END SUBROUTINE ROTATEV
!================================================================================================
SUBROUTINE FROMCART(R,LAT,LONG,POS)
!
!-----------------------------------------------------------------------
! CONVERT CARTESIAN COORDINATES POS(3)
! TO SPHERICAL COORDINATES R, LATITUDE, AND LONGITUDE (DEGREES)
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!------------------------------Arguments--------------------------------
!
REAL R, LAT, LONG, POS(3)
!
!---------------------------Local variables-----------------------------
!
real pi
!
!-----------------------------------------------------------------------
!
! pi=2.*ASIN(1.)
pi=3.141592653
R=SQRT(POS(1)*POS(1) + POS(2)*POS(2) + POS(3)*POS(3))
IF(R.EQ.0.)THEN
LAT=0.
LONG=0.
ELSE
LAT=ASIN(POS(3)*pi/180./R)
LONG=ATAN2(POS(2),POS(1))
LONG=LONG*180./pi
ENDIF
RETURN
END SUBROUTINE FROMCART
!================================================================================================
SUBROUTINE TOCART(R,LAT,LONG,POS)
!
!-----------------------------------------------------------------------
! CONVERT SPHERICAL COORDINATES R, LATITUDE, AND LONGITUDE (DEGREES)
! TO CARTESIAN COORDINATES POS(3)
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!------------------------------Arguments--------------------------------
!
REAL R, LAT, LONG, POS(3)
!
!---------------------------Local variables-----------------------------
!
real pi, stc, ctc, sf, cf
!
!-----------------------------------------------------------------------
!
! pi=2.*ASIN(1.)
pi=3.141592653
STC = SIN(LAT*pi/180.)
CTC = COS(LAT*pi/180.)
SF = SIN(LONG*pi/180.)
CF = COS(LONG*pi/180.)
POS(1) = R*CTC*CF
POS(2) = R*CTC*SF
POS(3) = R*STC
RETURN
END SUBROUTINE TOCART
!================================================================================================
SUBROUTINE ADJUST(ANGLE)
!
!-----------------------------------------------------------------------
! ADJUST AN ANGLE IN DEGREES TO BE IN RANGE OF 0 TO 360.
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!------------------------------Arguments--------------------------------
!
real angle
!
!-----------------------------------------------------------------------
!
10 CONTINUE
IF(ANGLE.LT.0.)THEN
ANGLE=ANGLE+360.
GOTO 10
ENDIF
20 CONTINUE
IF(ANGLE.GE.360.)THEN
ANGLE=ANGLE-360.
GOTO 20
ENDIF
RETURN
END SUBROUTINE ADJUST
!================================================================================================
INTEGER FUNCTION JULDAY_WAM(MM,ID,IYYY)
!
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
!
!------------------------------Arguments--------------------------------
!
integer mm, id, iyyy
!
!-----------------------------Parameters------------------------------
!
integer igreg
PARAMETER (IGREG=15+31*(10+12*1582))
!
!---------------------------Local variables-----------------------------
!
integer ja, jm, jy
!
!-----------------------------------------------------------------------
!
!!!compiler warning IF (IYYY.EQ.0) PAUSE 'There is no Year Zero.'
IF (IYYY.EQ.0) STOP 'There is no Year Zero.'
IF (IYYY.LT.0) IYYY=IYYY+1
IF (MM.GT.2) THEN
JY=IYYY
JM=MM+1
ELSE
JY=IYYY-1
JM=MM+13
ENDIF
JULDAY_WAM=INT(365.25*JY)+INT(30.6001*JM)+ID+1720995
IF (ID+31*(MM+12*IYYY).GE.IGREG) THEN
JA=INT(0.01*JY)
JULDAY_WAM=JULDAY_WAM+2-JA+INT(0.25*JA)
ENDIF
RETURN
END FUNCTION JULDAY_WAM
!================================================================================================
FUNCTION MLT(MagLong)
!
!-----------------------------------------------------------------------
! given magnetic longitude in degrees, return Magnetic Local Time
! assuming that TRANS has been called with the date & time to calculate
! the rotation matrices.
!
! btf 11/06/03:
! Call sub adjust instead of referencing it as a function
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only: CX=>CX,ST=>ST,CT=>CT,AM=>AM
implicit none
!
!-----------------------------Return Value------------------------------
!
real mlt
!
!-------------------------------Commons---------------------------------
!
! real cx, st, ct, am
! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
!
!------------------------------Arguments--------------------------------
!
REAL MagLong
!
!---------------------------Local variables-----------------------------
!
REAL angle, rotangle
!
!-----------------------------------------------------------------------
!
RotAngle=CX(7)
! MLT=ADJUST(Maglong+RotAngle+180.)/15.
angle = Maglong+RotAngle+180.
call adjust(angle)
mlt = angle/15.
RETURN
END FUNCTION MLT
!================================================================================================
FUNCTION MagLong(MLT)
!
!-----------------------------------------------------------------------
! return magnetic longitude in degrees, given Magnetic Local Time
! assuming that TRANS has been called with the date & time to calculate
! the rotation matrices.
!
! btf 11/06/03:
! Call sub adjust instead of referencing it as a function
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only:CX=>CX,ST=>ST,CT=>CT,AM=>AM
implicit none
!
!-----------------------------Return Value------------------------------
!
real MagLong
!
!-------------------------------Commons---------------------------------
!
! real cx, st, ct, am
! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
!
!------------------------------Arguments--------------------------------
!
REAL MLT
!
!---------------------------Local variables-----------------------------
!
REAL angle, rotangle
!
!-----------------------------------------------------------------------
!
RotAngle=CX(7)
angle=MLT*15.-RotAngle-180.
! MagLong=ADJUST(angle)
call adjust(angle)
MagLong = angle
RETURN
END FUNCTION MagLong
!================================================================================================
SUBROUTINE SunLoc(SunLat,SunLong)
!
!-----------------------------------------------------------------------
! Return latitude and longitude of sub-solar point.
! Assumes that TRANS has previously been called with the
! date & time to calculate the rotation matrices.
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only:CX=>CX,ST=>ST,CT=>CT,AM=>AM
implicit none
!
!-------------------------------Commons---------------------------------
!
! real cx, st, ct, am
! COMMON/TRANSDAT/CX(9),ST(6),CT(6),AM(3,3,11)
!
!------------------------------Arguments--------------------------------
!
Real SunLat,SunLong
!
!-----------------------------------------------------------------------
!
SunLong=CX(9)
SunLat=CX(8)
RETURN
END SUBROUTINE SunLoc
!================================================================================================
SUBROUTINE GECMP (AMLA,RMLT,ET,EP)
!
!-----------------------------------------------------------------------
! Get Electric field components for the Weimer electrostatic
! potential model. Before use, first load coefficients (CALL
! READCOEF) and initialize model conditions (CALL SETMODEL).
!
! INPUTS:
! AMLA = Absolute value of magnetic latitude (deg)
! RMLT = Magnetic local time (hours).
! RETURNS:
! ET = Etheta (magnetic equatorward*) E field component (V/m)
! EP = Ephi (magnetic eastward) E field component (V/m)
!
! * ET direction is along the magnetic meridian away from the
! current hemisphere; i.e., when ET > 0, the direction is
! southward when RMLA > 0
! northward when RMLA < 0
!
! NCAR addition (Jan 97). R.Barnes
!-----------------------------------------------------------------------
!
! use shr_kind_mod, only: r8 => shr_kind_r8
use efield, only: ALAMN =>ALAMN,ALAMX=>ALAMX,ALAMR=>ALAMR,
&STPD=>STPD,STP2=>STP2,CSTP=>CSTP,SSTP=>SSTP
implicit none
!
!-------------------------------Commons---------------------------------
!
! CECMP contains constants initialized in READCOEF
! real alamn, alamx, alamr, stpd, stp2, cstp, sstp
! COMMON /CECMP/ ALAMN,ALAMX,ALAMR,STPD,STP2,CSTP,SSTP
!
!------------------------------Arguments--------------------------------
!
real amla, rmlt, et, ep
!
!-----------------------------Parameters------------------------------
!
real d2r, r2d
PARAMETER ( D2R = 0.0174532925199432957692369076847 ,
& R2D = 57.2957795130823208767981548147)
!
!---------------------------Local variables-----------------------------
!
real p1, p2
real xmlt, xmlt1, kpol, dphi, amla1
!
!-------------------------External Functions----------------------------
!
real epotval
external epotval
!
!-----------------------------------------------------------------------
!
ET = -99999.
EP = -99999.
IF (AMLA .LT. 0.) GO TO 100
! Calculate -(latitude gradient) by stepping 10 km along the
! meridian in each direction (flipping coordinates when going
! over pole to keep lat <= 90).
KPOL = 0
XMLT = RMLT
10 XMLT1 = XMLT
AMLA1 = AMLA + STPD
IF (AMLA1 .GT. 90.) THEN
AMLA1 = 180. - AMLA1
XMLT1 = XMLT1 + 12.
ENDIF
P1 = EPOTVAL (AMLA1 ,XMLT1)
P2 = EPOTVAL (AMLA-STPD,XMLT )
IF (KPOL .EQ. 1) GO TO 20
ET = (P1 - P2) / STP2
! Calculate -(lon gradient). For most latitudes, step along a
! great circle. However, limit minimum latitude to the model
! minimum (distorting the path onto a latitude line). Also,
! avoid a divide by zero at the pole avoid by using Art's trick
! where Ephi(90,lon) = Etheta(90,lon+90)
IF (AMLA .LT. ALAMX) THEN
AMLA1 = MAX (ASIN(SIN(AMLA*D2R)*CSTP) , ALAMR)
DPHI = ASIN (SSTP/SIN(AMLA1))*R2D
AMLA1 = AMLA1*R2D
P1 = EPOTVAL (AMLA1,XMLT+DPHI)
P2 = EPOTVAL (AMLA1,XMLT-DPHI)
ELSE
AMLA = 90.
XMLT = XMLT + 6.
KPOL = 1
GO TO 10
ENDIF
20 EP = (P2 - P1) / STP2
IF (KPOL .EQ. 1) EP = -EP
! Below model minimum lat, the potential is value at min lat
IF (AMLA .LT. ALAMN) THEN
ET = 0.
EP = EP * COS(ALAMR)/COS(AMLA*D2R)
ENDIF
100 RETURN
END SUBROUTINE GECMP
!=====================================================================
subroutine svdcmp( a, m, n, mp, np, w, v )
!-------------------------------------------------------------------------
! purpose: singular value decomposition
!
! method:
! given a matrix a(1:m,1:n), with physical dimensions mp by np,
! this routine computes its singular value decomposition,
! the matrix u replaces a on output. the
! diagonal matrix of singular values w is output as a vector
! w(1:n). the matrix v (not the transpose v^t) is output as
! v(1:n,1:n).
!
! author: a. maute dec 2003
! (* copyright (c) 1985 numerical recipes software -- svdcmp *!
! from numerical recipes 1986 pp. 60 or can be find on web-sites
!-------------------------------------------------------------------------
implicit none
integer, parameter :: nmax = 1600
!-------------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------------
integer, intent(in) :: m
integer, intent(in) :: n
integer, intent(in) :: mp
integer, intent(in) :: np
real, intent(inout) :: a(mp,np)
real, intent(out) :: v(np,np)
real, intent(out) :: w(np)
!-------------------------------------------------------------------------
! ... local variables
!-------------------------------------------------------------------------
integer :: i, its, j, k, l, nm
real :: anorm
real :: c
real :: f
real :: g
real :: h
real :: s
real :: scale
real :: x, y, z
real :: rv1(nmax)
logical :: cnd1
logical :: cnd2
g = 0.0
scale = 0.0
anorm = 0.0
do i = 1,n !loop1
l = i + 1
rv1(i) = scale*g
g = 0.0
s = 0.0
scale = 0.0
if( i <= m ) then
do k = i,m
scale = scale + abs(a(k,i))
end do
if( scale /= 0.0 ) then
do k = i,m
a(k,i) = a(k,i)/scale
s = s + a(k,i)*a(k,i)
end do
f = a(i,i)
g = -sign(sqrt(s),f)
h = f*g - s
a(i,i) = f - g
if( i /= n ) then
do j = l,n
s = 0.0
do k = i,m
s = s + a(k,i)*a(k,j)
end do
f = s/h
do k = i,m
a(k,j) = a(k,j) + f*a(k,i)
end do
end do
end if
do k = i,m
a(k,i) = scale*a(k,i)
end do
endif
endif
w(i) = scale *g
g = 0.0
s = 0.0
scale = 0.0
if( i <= m .and. i /= n ) then
do k = l,n
scale = scale + abs(a(i,k))
end do
if( scale /= 0.0 ) then
do k = l,n
a(i,k) = a(i,k)/scale
s = s + a(i,k)*a(i,k)
end do
f = a(i,l)
g = -sign(sqrt(s),f)
h = f*g - s
a(i,l) = f - g
do k = l,n
rv1(k) = a(i,k)/h
end do
if( i /= m ) then
do j = l,m
s = 0.0
do k = l,n
s = s + a(j,k)*a(i,k)
end do
do k = l,n
a(j,k) = a(j,k) + s*rv1(k)
end do
end do
end if
do k = l,n
a(i,k) = scale*a(i,k)
end do
end if
end if
anorm = max( anorm,(abs(w(i)) + abs(rv1(i))) )
enddo !loop1
do i = n,1,-1
if( i < n ) then
if( g /= 0.0 ) then
do j = l,n
v(j,i) = (a(i,j)/a(i,l))/g
end do
do j = l,n
s = 0.0
do k = l,n
s = s + a(i,k)*v(k,j)
end do
do k = l,n
v(k,j) = v(k,j) + s*v(k,i)
end do
end do
end if
do j = l,n
v(i,j) = 0.0
v(j,i) = 0.0
end do
end if
v(i,i) = 1.0
g = rv1(i)
l = i
end do
do i = n,1,-1
l = i + 1
g = w(i)
if( i < n ) then
do j = l,n
a(i,j) = 0.0
end do
end if
if( g /= 0.0 ) then
g = 1./g
if( i /= n ) then
do j = l,n
s = 0.0
do k = l,m
s = s + a(k,i)*a(k,j)
end do
f = (s/a(i,i))*g
do k = i,m
a(k,j) = a(k,j) + f*a(k,i)
end do
end do
end if
do j = i,m
a(j,i) = a(j,i)*g
end do
else
do j = i,m
a(j,i) = 0.0
end do
end if
a(i,i) = a(i,i) + 1.0
end do
do k = n,1,-1
do its = 1,30 !loop2
do l = k,1,-1
nm = l - 1
cnd1 = abs( rv1(l) ) + anorm == anorm
if( cnd1 ) then
cnd2 = .false.
exit
end if
cnd2 = abs( w(nm) ) + anorm == anorm
if( cnd2 ) then
cnd1 = .true.
exit
else if( l == 1 ) then
cnd1 = .true.
cnd2 = .true.
end if
end do
if( cnd2 ) then
c = 0.0
s = 1.0
do i = l,k
f = s*rv1(i)
if( (abs(f) + anorm) /= anorm ) then
g = w(i)
h = sqrt(f*f + g*g)
w(i) = h
h = 1.0/h
c = (g*h)
s = -(f*h)
do j = 1,m
y = a(j,nm)
z = a(j,i)
a(j,nm) = (y*c) + (z*s)
a(j,i) = -(y*s) + (z*c)
end do
end if
end do
end if
if( cnd1 ) then
z = w(k)
if( l == k ) then
if( z < 0.0 ) then
w(k) = -z
do j = 1,n
v(j,k) = -v(j,k)
end do
end if
! exit loop2
go to 20
end if
end if
x = w(l)
nm = k - 1
y = w(nm)
g = rv1(nm)
h = rv1(k)
f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y)
g = sqrt( f*f + 1.0 )
f = ((x - z)*(x + z) + h*((y/(f + sign(g,f))) - h))/x
c = 1.0
s = 1.0
do j = l,nm
i = j + 1
g = rv1(i)
y = w(i)
h = s*g
g = c*g
z = sqrt( f*f + h*h )
rv1(j) = z
c = f/z
s = h/z
f = (x*c)+(g*s)
g = -(x*s)+(g*c)
h = y*s
y = y*c
do nm = 1,n
x = v(nm,j)
z = v(nm,i)
v(nm,j) = (x*c)+(z*s)
v(nm,i) = -(x*s)+(z*c)
end do
z = sqrt( f*f + h*h )
w(j) = z
if( z /= 0.0 ) then
z = 1.0/z
c = f*z
s = h*z
end if
f = (c*g)+(s*y)
x = -(s*g)+(c*y)
do nm = 1,m
y = a(nm,j)
z = a(nm,i)
a(nm,j) = (y*c)+(z*s)
a(nm,i) = -(y*s)+(z*c)
end do
end do
rv1(l) = 0.0
rv1(k) = f
w(k) = x
end do !loop2
20 continue
end do
end subroutine svdcmp
!-------------------------------------------------------------------------
! purpose: solves a*x = b
!
! method:
! solves a*x = b for a vector x, where a is specified by the arrays
! u,w,v as returned by svdcmp. m and n
! are the logical dimensions of a, and will be equal for square matrices.
! mp and np are the physical dimensions of a. b(1:m) is the input right-hand
! side. x(1:n) is the output solution vector. no input quantities are
! destroyed, so the routine may be called sequentially with different b
!
! author: a. maute dec 2002
! (* copyright (c) 1985 numerical recipes software -- svbksb *!
! from numerical recipes 1986 pp. 57 or can be find on web-sites
!-------------------------------------------------------------------------
subroutine svbksb( u, w, v, m, n, mp, np, b, x )
!-------------------------------------------------------------------------
! ... dummy arguments
!-------------------------------------------------------------------------
implicit none
integer, parameter :: nmax = 1600
integer, intent(in) :: m
integer, intent(in) :: n
integer, intent(in) :: mp
integer, intent(in) :: np
real , intent(in) :: u(mp,np)
real , intent(in) :: w(np)
real , intent(in) :: v(np,np)
real , intent(in) :: b(mp)
real , intent(out) :: x(np)
!-------------------------------------------------------------------------
! ... local variables
!-------------------------------------------------------------------------
integer :: i, j, jj
real :: s
real :: tmp(nmax)
do j = 1,n
s = 0.
if( w(j) /= 0. ) then
do i = 1,m
s = s + u(i,j)*b(i)
end do
s = s/w(j)
endif
tmp(j) = s
end do
do j = 1,n
s = 0.
do jj = 1,n
s = s + v(j,jj)*tmp(jj)
end do
x(j) = s
end do
end subroutine svbksb