Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! Heelis empirical model to calculate high-latitude electric potential
3 : !-----------------------------------------------------------------------
4 : module heelis_mod
5 :
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use aurora_params, only: offc, dskofc, phin, phid, theta0
8 : use aurora_params, only: ctpoten, hpower, plevel, dskofa, rrad
9 :
10 : implicit none
11 : private
12 :
13 : public :: heelis_update
14 : public :: heelis_flwv32
15 :
16 : ! private data
17 : ! Auroral parameters (taken from aurora.F of timegcm):
18 : ! (dimension of 2 is for south,north hemispheres)
19 : !
20 : real(r8) :: &
21 : psim(2) ,& ! night convection entrance in MLT converted to radians (f(By))
22 : psie(2) ,& !
23 : pcen(2) ,& !
24 : phidp0(2) ,& !
25 : phidm0(2) ,& !
26 : phinp0(2) ,& !
27 : phinm0(2) ,& !
28 : rr1(2) !
29 :
30 : contains
31 : !-----------------------------------------------------------------------
32 8064 : subroutine heelis_update(max_ctpoten)
33 : use physconst, only: pi
34 : use mag_parms, only: get_mag_parms
35 :
36 : real(r8), optional, intent(in) :: max_ctpoten ! max cross cap potential kV
37 : !
38 : real(r8) :: rcp, rhp, arad
39 : real(r8) :: byloc ! local By; now is just a hook, and set to 0.
40 : integer, parameter :: isouth = 1
41 : integer, parameter :: inorth = 2
42 : real(r8), parameter :: h2deg = 15._r8 ! hour to degree
43 : real(r8), parameter :: dtr = pi/180._r8
44 :
45 : !
46 : ! This is called at every timestep because ctpoten may change with time.
47 : !
48 8064 : call get_mag_parms( hpower = hpower, ctpoten = ctpoten )
49 8064 : if (present(max_ctpoten)) then
50 0 : ctpoten = min(ctpoten,max_ctpoten)
51 : end if
52 :
53 : byloc = 0._r8
54 :
55 24192 : offc(:) = 1._r8*dtr
56 8064 : dskofc(:) = 0._r8
57 24192 : phin(:) = 180._r8*dtr
58 :
59 8064 : phid(isouth) = (9.39_r8 + 0.21_r8*byloc - 12._r8) * h2deg * dtr
60 8064 : phid(inorth) = (9.39_r8 - 0.21_r8*byloc - 12._r8) * h2deg * dtr
61 8064 : phin(isouth) = (23.50_r8 + 0.15_r8*byloc - 12._r8) * h2deg * dtr
62 8064 : phin(inorth) = (23.50_r8 - 0.15_r8*byloc - 12._r8) * h2deg * dtr
63 24192 : psim(:) = 0.44_r8 * ctpoten * 1000._r8
64 24192 : psie(:) = -0.56_r8 * ctpoten * 1000._r8
65 8064 : pcen(isouth) = (-0.168_r8 + 0.027_r8*byloc) * ctpoten * 1000._r8
66 8064 : pcen(inorth) = (-0.168_r8 - 0.027_r8*byloc) * ctpoten * 1000._r8
67 :
68 24192 : phidp0(:) = 90._r8*dtr
69 24192 : phidm0(:) = 90._r8*dtr
70 24192 : phinp0(:) = 90._r8*dtr
71 24192 : phinm0(:) = 90._r8*dtr
72 24192 : rr1(:) = -2.6_r8
73 24192 : theta0(:) = (-3.80_r8+8.48_r8*(ctpoten**0.1875_r8))*dtr
74 :
75 8064 : dskofa(:) = 0._r8
76 :
77 8064 : if( hpower >= 1.0_r8 ) then
78 8064 : plevel = 2.09_r8*log( hpower )
79 : else
80 0 : plevel = 0._r8
81 : end if
82 8064 : rhp = 14.20_r8 + 0.96_r8*plevel
83 8064 : rcp = -0.43_r8 + 9.69_r8 * (ctpoten**.1875_r8)
84 8064 : arad = max( rcp,rhp )
85 24192 : rrad(:) = arad*dtr
86 :
87 8064 : end subroutine heelis_update
88 :
89 : !-----------------------------------------------------------------------
90 387072 : subroutine heelis_flwv32(dlat,dlon,ratio,pi_in,iflag,nmlon,poten)
91 : !
92 : ! Calculate heelis potential at current magnetic latitude mlat.
93 : !
94 : !
95 : ! Args:
96 : real(r8), intent(in) :: pi_in
97 : integer, intent(in) :: nmlon
98 : integer, intent(inout) :: iflag(nmlon)
99 : real(r8),dimension(nmlon), intent(in) :: dlat,dlon,ratio
100 : real(r8),dimension(nmlon+1), intent(out) :: poten
101 : !
102 : ! Local:
103 : integer :: i,n,ihem
104 : real(r8),parameter :: eps=1.e-6_r8
105 : real(r8) :: &
106 : pi2,pih,sinthr1,psi(8),phirc,sinth0, &
107 : ofdc,cosofc(2),sinofc(2),aslonc(2), &
108 : phdpmx(2),phnpmx(2),phnmmx(2),phdmmx(2)
109 774144 : real(r8),dimension(nmlon) :: sinlat,coslat,sinlon,coslon,alon, &
110 774144 : colat,wk1,wk2,wk3,phifun,phifn2
111 774144 : integer :: ifn(nmlon)
112 774144 : real(r8) :: phi(nmlon,8)
113 : !
114 387072 : pi2 = 2.0_r8*pi_in
115 387072 : pih = 0.5_r8*pi_in
116 :
117 1161216 : do n=1,2
118 : !
119 774144 : ofdc = sqrt(offc(n)**2+dskofc(n)**2)
120 774144 : cosofc(n) = cos(ofdc)
121 774144 : sinofc(n) = sin(ofdc)
122 774144 : aslonc(n) = asin(dskofc(n)/ofdc)
123 : !
124 774144 : if (phin(n) < phid(n)) phin(n) = phin(n)+pi2 ! modifies aurora phin
125 774144 : phdpmx(n) = .5_r8*min(pi_in,(phin(n)-phid(n)))
126 774144 : phnpmx(n) = .5_r8*min(pi_in,(phid(n)-phin(n)+pi2))
127 774144 : phnmmx(n) = phdpmx(n)
128 1161216 : phdmmx(n) = phnpmx(n)
129 : enddo ! n=1,2
130 : !
131 : ! Set ihem=1,2 for South,North hemisphere:
132 : !
133 387072 : ihem = int( dlat(max0(1,nmlon/2))*2._r8/3.1416_r8 + 2._r8 )
134 387072 : sinth0 = sin(theta0(ihem))
135 : !
136 : ! Average amie results show r1=-2.6 for 11.3 degrees
137 : ! (0.1972 rad) beyond theta0.
138 : !
139 387072 : sinthr1 = sin(theta0(ihem)+0.1972_r8)
140 387072 : psi(1) = psie(ihem)
141 387072 : psi(3) = psim(ihem)
142 387072 : do n=2,4,2
143 774144 : psi(n) = psi(n-1)
144 : enddo ! n=2,4,2
145 1935360 : do n=1,4
146 1935360 : psi(n+4) = psi(n)
147 : enddo ! n=1,4
148 : !
149 : ! Transform to auroral circle coordinates:
150 : !
151 31352832 : do i=1,nmlon
152 30965760 : sinlat(i) = sin(abs(dlat(i)))
153 30965760 : coslat(i) = cos(dlat(i))
154 30965760 : sinlon(i) = sin(dlon(i)+aslonc(ihem))
155 30965760 : coslon(i) = cos(dlon(i)+aslonc(ihem))
156 : colat(i) = cosofc(ihem)*sinlat(i)-sinofc(ihem)*coslat(i)* &
157 30965760 : coslon(i)
158 : alon(i) = mod(atan2(sinlon(i)*coslat(i),sinlat(i)* &
159 : sinofc(ihem)+cosofc(ihem)*coslat(i)*coslon(i))- &
160 30965760 : aslonc(ihem)+3._r8*pi_in,pi2)-pi_in
161 30965760 : colat(i) = acos(colat(i))*sqrt(ratio(i))
162 : !
163 : ! Boundaries for longitudinal function:
164 : !
165 30965760 : wk1(i) = ((colat(i)-theta0(ihem))/theta0(ihem))**2
166 : phi(i,4)=phid(ihem)+eps-min(phidm0(ihem)+wk1(i)* &
167 30965760 : (pih-phidm0(ihem)),phdmmx(ihem))
168 : phi(i,5)=phid(ihem)-eps+min(phidp0(ihem)+wk1(I)* &
169 30965760 : (pih-phidp0(ihem)),phdpmx(ihem))
170 : phi(i,6)=phin(ihem)+eps-min(phinm0(ihem)+wk1(i)* &
171 30965760 : (pih-phinm0(ihem)),phnmmx(ihem))
172 : phi(i,7)=phin(ihem)-eps+min(phinp0(ihem)+wk1(i)* &
173 30965760 : (pih-phinp0(ihem)),phnpmx(ihem))
174 30965760 : phi(i,1)=phi(i,5)-pi2
175 30965760 : phi(i,2)=phi(i,6)-pi2
176 30965760 : phi(i,3)=phi(i,7)-pi2
177 30965760 : phi(i,8)=phi(i,4)+pi2
178 30965760 : phifun(i)=0._r8
179 30965760 : phifn2(i) = 0._r8
180 30965760 : if (colat(i)-theta0(ihem) >= 0._r8) then
181 20922624 : ifn(i) = 3
182 : else
183 10043136 : ifn(i) = 2
184 : endif
185 30965760 : if (iflag(i) == 1) iflag(i) = ifn(i)
186 : !
187 : ! Add ring current rotation to potential (phirc)
188 : !
189 30965760 : phirc = 0._r8
190 30965760 : wk2(i) = mod(alon(i)+phirc+2._r8*pi2+pi_in,pi2)-pi_in
191 31352832 : wk3(i) = mod(alon(i)+phirc+3._r8*pi2,pi2)-pi_in
192 : enddo ! i=1,nmlon
193 : !
194 : ! Longitudinal variation:
195 : !
196 3096576 : do n=1,7
197 219856896 : do i=1,nmlon
198 650280960 : phifun(i)=phifun(i)+.25_r8*(psi(n)+psi(n+1)+(psi(n)- &
199 : psi(n+1))*cos(mod(pi_in*(wk2(i)-phi(i,n))/(phi(i,n+1)- &
200 : phi(i,n)),pi2)))*(1._r8-sign(1._r8,(wk2(i)-phi(i,n))* &
201 650280960 : (wk2(i)-phi(i,n+1))))
202 : phifn2(i)=phifn2(i)+.25_r8*(psi(n)+psi(n+1)+(psi(n)- &
203 : psi(n+1))*cos(mod(pi_in*(wk3(i)-phi(i,n))/(phi(i,n+1)- &
204 : phi(i,n)),pi2)))*(1._r8-sign(1._r8,(wk3(i)-phi(i,n))* &
205 219469824 : (wk3(i)-phi(i,n+1))))
206 : enddo
207 : enddo
208 : !
209 : ! Evaluate total potential:
210 : !
211 :
212 31352832 : do i=1,nmlon
213 31352832 : if (iflag(i)==2) then
214 10043136 : poten(i) = (2._r8*(pcen(ihem)-phifun(i))+(phifun(i)-phifn2(i))* &
215 : 0.75_r8)*(colat(i)/theta0(ihem))**3 + &
216 : (1.5_r8*(phifun(i)+phifn2(i))-3._r8*pcen(ihem))*(colat(i)/ &
217 : theta0(ihem))**2 + 0.75_r8*(phifun(i)-phifn2(i))*(colat(i)/ &
218 10043136 : theta0(ihem)) + pcen(ihem)
219 : else
220 20922624 : poten(i) = phifun(i)*(max(sin(colat(i)),sinth0)/sinth0)**rr1(ihem)* &
221 20922624 : exp(7._r8*(1._r8-max(sin(colat(i)),sinthr1)/sinthr1))
222 : endif
223 : enddo
224 :
225 8064 : end subroutine heelis_flwv32
226 : !-----------------------------------------------------------------------
227 : end module heelis_mod
|
