Line data Source code
1 :
2 : module efield
3 : !------------------------------------------------------------------------------
4 : ! description: calculates the electric potential for a given year,
5 : ! day of year, UT, F10.7, K_p
6 : ! - low/midlatitudes electric potential is from an empirical model from
7 : ! L.Scherliess ludger@gaim.cass.usu.edu
8 : ! - high latitude electric potential is Heelis
9 : ! - the transition zone is smoothed
10 : ! - output is the horizontal global electric field in magnetic coordinates direction
11 : ! at every magnetic local time grid point expressed in degrees (0 deg-0MLT; 360deg 24 MLT)
12 : !
13 : ! input
14 : ! integer :: iday, ! day number of year
15 : ! iyear ! year
16 : ! real(r8):: ut, ! universal time
17 : ! F10.7, ! solar flux (see ionosphere module)
18 : ! output
19 : ! real(r8) :: &
20 : ! ed1(0:nmlon,0:nmlat), & ! zonal electric field Ed1 [V/m]
21 : ! ed2(0:nmlon,0:nmlat) ! meridional electric field Ed2/sin I_m [V/m]
22 : !
23 : ! notes:
24 : !
25 : ! - !to be done (commented out): input S_a F10.7/ Kp from WACCM and calculate B_z
26 : ! from these inputs
27 : ! - assume regular geomagnetic grid
28 : ! - uses average year 365.24 days/year 30.6001 day/mo s. Weimer
29 : ! - get_tilt works only for iyear >= 1900
30 : ! - fixed parameters: B_z, B_y units nT CHANGE THIS
31 : ! F10.7
32 : ! - we assume that the reference height is 300km for the emperical potential model
33 : ! - as a first approximation the electric field is constant in height
34 : ! WATCH what is the upper boundary condition in WACCM
35 : ! - for all the calculation done here we set the reference height to the same
36 : ! value as in tiegcm (hr=130km)
37 : ! - 12/15/03 input value iseasav : replaced by day -> month and day of month
38 : ! - 12/15/03 S_aM calculated according to Scherliess draft paper and added
39 : ! S_aM(corrected) = 90*(S_aM+1) to get variation in fig 1 Scherliess draft
40 : !
41 : ! Author: A. Maute Dec 2003 am 12/30/03
42 : !------------------------------------------------------------------------------
43 :
44 : use shr_kind_mod, only: r8 => shr_kind_r8
45 : use physconst, only: pi
46 : use cam_abortutils, only: endrun
47 : use cam_logfile, only: iulog
48 : use heelis_mod, only: heelis_flwv32, heelis_update
49 : use solar_parms_data, only: f107d=>solar_parms_f107
50 : use spmd_utils, only: masterproc
51 : use infnan, only: nan, assignment(=)
52 :
53 : implicit none
54 :
55 : public :: efield_init, & ! interface routine
56 : get_efield ! interface routine
57 : public :: ed1, & ! zonal electric field Ed1 [V/m]
58 : ed2, & ! meridional electric field Ed2 [V/m]
59 : potent, & ! electric potential [V]
60 : nmlon, nmlat, & ! dimension of mag. grid
61 : dlatm, dlonm, & ! grid spacing of mag. grid
62 : ylonm, ylatm ! magnetic longitudes/latitudes (deg)
63 :
64 : private
65 :
66 : integer :: &
67 : iday, & ! day number of year
68 : iyear, & ! year
69 : iday_m, & ! day of month
70 : imo ! month
71 : real(r8) :: ut ! universal time
72 :
73 : !------------------------------------------------------------------------------
74 : ! mag. grid dimensions (assumed resolution of 2deg)
75 : !------------------------------------------------------------------------------
76 : integer, parameter :: &
77 : nmlon = 180, & ! mlon
78 : nmlat = 90, & ! mlat
79 : nmlath= nmlat/2 ! mlat/2
80 :
81 : integer, parameter :: &
82 : nmlon1f = nmlon/4, & ! 1 fourth mlon
83 : nmlon2f = nmlon/2, & ! 2 fourths mlon
84 : nmlon3f = 3*nmlon/4 ! 3 fourths mlon
85 :
86 : real(r8) :: &
87 : ylatm(0:nmlat), & ! magnetic latitudes (deg)
88 : ylonm(0:nmlon), & ! magnetic longitudes (deg)
89 : dlonm, & ! delon lon grid spacing
90 : dlatm ! delat lat grid spacing
91 :
92 : !------------------------------------------------------------------------------
93 : ! array on magnetic grid:
94 : !------------------------------------------------------------------------------
95 : real(r8), protected :: &
96 : potent(0:nmlon,0:nmlat), & ! electric potential [V]
97 : ed1(0:nmlon,0:nmlat), & ! zonal electric field Ed1 [V/m]
98 : ed2(0:nmlon,0:nmlat) ! meridional electric field Ed2/sin I_m [V/m]
99 :
100 : real(r8) :: &
101 : day ! iday+ut
102 :
103 : logical, parameter :: iutav=.false. ! .true. means UT-averaging
104 : ! .false. means no UT-averaging
105 : !------------------------------------------------------------------------------
106 : ! boundary for high-lat potential
107 : !------------------------------------------------------------------------------
108 : real(r8), parameter :: bnd_hgh = 44._r8 ! colat. [deg]
109 : integer :: nmlat_hgh
110 :
111 : !------------------------------------------------------------------------------
112 : ! flag for choosing factors for empirical low latitude model
113 : !------------------------------------------------------------------------------
114 : integer, parameter :: iseasav = 0 ! flag for season
115 :
116 : !------------------------------------------------------------------------------
117 : ! constants:
118 : !------------------------------------------------------------------------------
119 : real(r8), parameter :: &
120 : r_e = 6.371e6_r8, & ! radius_earth [m] (same as for apex.F90)
121 : h_r = 130.0e3_r8, & ! reference height [m] (same as for apex.F90)
122 : dy2yr= 365.24_r8 ! day per avg. year
123 :
124 : real(r8) :: &
125 : rtd , & ! radians -> deg
126 : dtr, & ! deg -> radians
127 : sqr2, &
128 : dy2rd, & ! 2*pi/365.24 average year
129 : sinIm_mag(0:nmlat) ! sinIm
130 :
131 : integer :: jmin, jmax ! latitude index for interpolation of
132 : ! northward e-field ed2 at mag. equator
133 :
134 : !------------------------------------------------------------------------------
135 : ! for spherical harmonics
136 : !------------------------------------------------------------------------------
137 : integer, parameter :: &
138 : nm = 19, &
139 : mm = 18, &
140 : nmp = nm + 1, &
141 : mmp = mm + 1
142 :
143 : real(r8) :: r(0:nm,0:mm) ! R_n^m
144 : real(r8) :: pmopmmo(0:mm) ! sqrt(1+1/2m)
145 :
146 : !------------------------------------------------------------------------------
147 : ! index for factors f_m(mlt),f_l(UT),f_-k(d)
148 : !------------------------------------------------------------------------------
149 : integer, parameter :: ni = 1091 ! for n=12 m=-18:18
150 : integer :: imax ! max number of index
151 : integer,dimension(0:ni) :: &
152 : kf, &
153 : lf, &
154 : mf, &
155 : nf, &
156 : jf
157 : real(r8) :: ft(1:3,0:2) ! used for f_-k(season,k)
158 :
159 : real(r8) :: a_klnm(0:ni) ! A_klm
160 : real(r8) :: a_lf(0:ni) ! A_klmn^lf for minimum &
161 : real(r8) :: a_hf(0:ni) ! A_klmn^hf for maximum
162 :
163 : !------------------------------------------------------------------------------
164 : ! high_latitude boundary
165 : !------------------------------------------------------------------------------
166 : real(r8), parameter :: &
167 : lat_sft = 54._r8 ! shift of highlat_bnd to 54 deg
168 : integer :: ilat_sft ! index of shift for high latitude boundary
169 : integer, parameter :: nmax_sin = 2 ! max. wave number to be represented
170 : logical, parameter :: debug =.false.
171 :
172 : real(r8) :: epotential_max = huge(1._r8) ! max cross cap potential kV
173 :
174 : integer :: ip1f(0:nmlon)=-huge(1)
175 : integer :: ip2f(0:nmlon)=-huge(1)
176 : integer :: ip3f(0:nmlon)=-huge(1)
177 :
178 : contains
179 :
180 1536 : subroutine efield_init(efield_lflux_file, efield_hflux_file, efield_potential_max)
181 : !-----------------------------------------------------------------------
182 : ! Purpose: read in and set up coefficients needed for electric field
183 : ! calculation (independent of time & geog. location)
184 : !
185 : ! Method:
186 : !
187 : ! Author: A. Maute Dec 2003 am 12/17/03
188 : !-----------------------------------------------------------------------
189 : character(len=*), intent(in) :: efield_lflux_file
190 : character(len=*), intent(in) :: efield_hflux_file
191 : real(r8), intent(in) :: efield_potential_max ! cross cap electric potential maximum
192 :
193 : integer :: i
194 : real(r8) :: nanval
195 1536 : nanval=nan
196 :
197 25440768 : potent = nanval
198 25440768 : ed1 = nanval
199 25440768 : ed2 = nanval
200 :
201 1536 : call constants ! calculate constants
202 : !-----------------------------------------------------------------------
203 : ! low/midlatitude potential from Scherliess model
204 : !-----------------------------------------------------------------------
205 1536 : call read_acoef (efield_lflux_file, efield_hflux_file) ! read in A_klnm for given S_aM
206 1536 : call index_quiet ! set up index for f_m(mlt),f_l(UT),f_-k(d)
207 1536 : call prep_fk ! set up the constant factors for f_k
208 1536 : call prep_pnm ! set up the constant factors for P_n^m & dP/d phi
209 :
210 1536 : epotential_max = efield_potential_max
211 :
212 279552 : do i = 0,nmlon
213 278016 : ip1f(i) = i + nmlon1f
214 278016 : if( ip1f(i) > nmlon ) then
215 69120 : ip1f(i) = ip1f(i) - nmlon
216 : end if
217 278016 : ip2f(i) = i + nmlon2f
218 278016 : if( ip2f(i) > nmlon ) then
219 138240 : ip2f(i) = ip2f(i) - nmlon
220 : end if
221 278016 : ip3f(i) = i + nmlon3f
222 279552 : if( ip3f(i) > nmlon ) then
223 207360 : ip3f(i) = ip3f(i) - nmlon
224 : end if
225 : end do
226 :
227 1536 : end subroutine efield_init
228 :
229 16128 : subroutine get_efield
230 : !-----------------------------------------------------------------------
231 : ! Purpose: calculates the global electric potential field on the
232 : ! geomagnetic grid (MLT in deg) and derives the electric field
233 : !
234 : ! Method:
235 : !
236 : ! Author: A. Maute Dec 2003 am 12/17/03
237 : !-----------------------------------------------------------------------
238 :
239 : use time_manager, only : get_curr_calday, get_curr_date
240 : use mo_apex, only : geomag_year
241 :
242 : integer :: tod ! time of day [s]
243 : character(len=80) :: errmsg
244 :
245 : !-----------------------------------------------------------------------
246 : ! get current calendar day of year & date components
247 : ! valid at end of current timestep
248 : !-----------------------------------------------------------------------
249 16128 : iday = get_curr_calday() ! day of year
250 16128 : call get_curr_date (iyear,imo,iday_m,tod) ! year, time of day [sec]
251 16128 : iyear = int(geomag_year)
252 :
253 16128 : if( iyear < 1900 ) then
254 0 : write(errmsg,"(/,'>>> get_efield: year < 1900 not possible: year=',i5)") iyear
255 0 : call endrun(errmsg)
256 : end if
257 :
258 16128 : ut = tod/3600._r8 ! UT of day [hour]
259 :
260 : !-----------------------------------------------------------------------
261 : ! ajust S_a
262 : !-----------------------------------------------------------------------
263 16128 : call adj_S_a
264 : !-----------------------------------------------------------------------
265 : ! calculate global electric potential
266 : !-----------------------------------------------------------------------
267 16128 : call GlobalElPotential
268 :
269 : !-----------------------------------------------------------------------
270 : ! calculate derivative of global electric potential
271 : !-----------------------------------------------------------------------
272 16128 : call DerivPotential
273 :
274 16128 : end subroutine get_efield
275 :
276 16128 : subroutine GlobalElPotential
277 : !-----------------------------------------------------------------------
278 : ! Purpose: calculates the global electric potential field on the
279 : ! geomagnetic grid (MLT in deg)
280 : !
281 : ! Method: rewritten code from Luedger Scherliess (11/20/99 LS)
282 : ! routine to calculate the global electric potential in magnetic
283 : ! Apex coordinates (Latitude and MLT).
284 : ! High Latitude Model is Heelis
285 : ! Midlatitude model is Scherliess 1999.
286 : ! Interpolation in a transition region at about 60 degree
287 : ! magnetic apex lat
288 : !
289 : ! Author: A. Maute Dec 2003 am 12/17/03
290 : !-----------------------------------------------------------------------
291 :
292 : !-----------------------------------------------------------------------
293 : ! local variables
294 : !-----------------------------------------------------------------------
295 : integer :: ilon, ilat, idlat
296 : integer :: ihlat_bnd(0:nmlon) ! high latitude boundary
297 : integer :: itrans_width(0:nmlon) ! width of transition zone
298 : real(r8) :: mlat, pot
299 : real(r8) :: pot_midlat(0:nmlon,0:nmlat) ! potential from L. Scherliess model
300 : real(r8) :: pot_highlat(0:nmlon,0:nmlat) ! potential from Heelis
301 : real(r8) :: pot_highlats(0:nmlon,0:nmlat) ! smoothed potential from Heelis
302 : real(r8) :: poten(nmlon+1)
303 : integer,dimension(nmlon) :: iflag
304 : real(r8),dimension(nmlon) :: dlat,dlon,ratio
305 :
306 : !-----------------------------------------------------------------------
307 : ! convert to date and day
308 : !-----------------------------------------------------------------------
309 16128 : day = iday + ut/24._r8
310 :
311 : !-----------------------------------------------------------------------
312 : ! low/mid-latitude electric potential - empirical model Scherliess 1999
313 : !-----------------------------------------------------------------------
314 : !$omp parallel do private(ilat, ilon, mlat, pot)
315 758016 : do ilat = 0,nmlath ! Calculate only for one magn. hemisphere
316 741888 : mlat = ylatm(ilat) ! mag. latitude
317 135039744 : do ilon = 0,nmlon ! lon. loop
318 134281728 : call efield_mid( mlat, ylonm(ilon), pot )
319 134281728 : pot_midlat(ilon,ilat+nmlath) = pot ! SH/NH symmetry
320 269305344 : pot_midlat(ilon,nmlath-ilat) = pot
321 : end do
322 : end do
323 :
324 : !-----------------------------------------------------------------------
325 : ! high-latitude potential from Heelis model
326 : !-----------------------------------------------------------------------
327 16128 : call heelis_update(max_ctpoten=epotential_max)
328 :
329 2919168 : ratio(:) = 1._r8
330 :
331 1483776 : do ilat = 0,nmlat
332 265644288 : iflag(:) = 1
333 265644288 : dlat(:) = (90._r8-ylatm(ilat))*dtr
334 265644288 : dlon(:) = (ylonm(1:nmlon)-180._r8)*dtr
335 :
336 1467648 : if (ylatm(ilat) > 44._r8) then
337 1096704 : call heelis_flwv32(dlat,dlon,ratio,pi,iflag,nmlon,poten(:))
338 : else
339 370944 : poten(:) = 0._r8
340 : endif
341 :
342 265644288 : pot_highlat(1:nmlon,ilat) = poten(1:nmlon) ! volts
343 265644288 : pot_highlat(1:nmlon,nmlat-ilat) = poten(1:nmlon)
344 265644288 : pot_highlats(1:nmlon,ilat) = poten(1:nmlon)
345 265660416 : pot_highlats(1:nmlon,nmlat-ilat) = poten(1:nmlon)
346 :
347 : end do
348 1483776 : pot_highlat(0,:) = pot_highlat(nmlon,:)
349 1483776 : pot_highlats(0,:) = pot_highlats(nmlon,:)
350 :
351 : !-----------------------------------------------------------------------
352 : ! weighted smoothing of high latitude potential
353 : !-----------------------------------------------------------------------
354 16128 : idlat = 2 ! smooth over -2:2 = 5 grid points
355 16128 : call pot_latsmo( pot_highlats, idlat )
356 : !-----------------------------------------------------------------------
357 : ! calculate the height latitude bounday ihl_bnd
358 : ! 1. calculate E field from high-lat potential model
359 : ! boundary is set where the total electric field exceeds
360 : ! 0.015 V/m (corresp. approx. to 300 m/s)
361 : ! 2. moved halfways to 54 deg
362 : ! output : index 0-pole nmlath-equator
363 : !-----------------------------------------------------------------------
364 16128 : potent = pot_highlat ! need efield components (ed1,ed2) to determine ihlat_bnd
365 16128 : call DerivPotential()
366 16128 : call highlat_getbnd( ihlat_bnd )
367 : !-----------------------------------------------------------------------
368 : ! 3. adjust high latitude boundary sinusoidally
369 : ! calculate width of transition zone
370 : !-----------------------------------------------------------------------
371 16128 : call bnd_sinus( ihlat_bnd, itrans_width )
372 : !-----------------------------------------------------------------------
373 : ! 4. ajust high latitude potential to low latitude potential
374 : !-----------------------------------------------------------------------
375 16128 : call highlat_adjust( pot_highlats, pot_highlat, pot_midlat, ihlat_bnd )
376 : !-----------------------------------------------------------------------
377 : ! interpolation of high and low/midlatitude potential in the
378 : ! transition zone and put it into global potent array
379 : !-----------------------------------------------------------------------
380 16128 : call interp_poten( pot_highlats, pot_highlat, pot_midlat, ihlat_bnd, itrans_width)
381 : !-----------------------------------------------------------------------
382 : ! potential weighted smoothing in latitude
383 : !-----------------------------------------------------------------------
384 16128 : idlat = 2 ! smooth over -2:2 = 5 grid points
385 16128 : call pot_latsmo2( potent, idlat )
386 : !-----------------------------------------------------------------------
387 : ! potential smoothing in longitude
388 : !-----------------------------------------------------------------------
389 16128 : idlat = nmlon/48 ! smooth over -idlat:idlat grid points
390 16128 : call pot_lonsmo( potent, idlat )
391 :
392 16128 : end subroutine GlobalElPotential
393 :
394 271482624 : subroutine ff( ph, mt, f )
395 : !-----------------------------------------------------------------------
396 : ! Purpose: calculate F for normalized associated Legendre polynomial P_n^m
397 : ! Ref.: Richmond J.Atm.Ter.Phys. 1974
398 : !
399 : ! Method: f_m(phi) = sqrt(2) sin(m phi) m > 0
400 : ! = 1 m = 0
401 : ! = sqrt(2) cos(m phi) m < 0
402 : !
403 : ! Author: A. Maute Nov 2003 am 11/18/03
404 : !-----------------------------------------------------------------------
405 :
406 : implicit none
407 :
408 : !-----------------------------------------------------------------------
409 : ! dummy arguments
410 : !-----------------------------------------------------------------------
411 : integer,intent(in) :: mt
412 : real(r8),intent(in) :: ph ! geo. longitude of 0SLT (ut*15)
413 : real(r8),intent(out) :: f(-mt:mt)
414 :
415 : !-----------------------------------------------------------------------
416 : ! local variables
417 : !-----------------------------------------------------------------------
418 : integer :: m, mmo
419 : real(r8) :: sp, cp
420 :
421 271482624 : sp = sin( ph/rtd )
422 271482624 : cp = cos( ph/rtd )
423 271482624 : f(0) = 1.e0_r8
424 :
425 271482624 : f(-1) = sqr2*cp
426 271482624 : f(1) = sqr2*sp
427 2691472896 : do m = 2,mt
428 2419990272 : mmo = m - 1
429 2419990272 : f(m) = f(-mmo)*sp + cp*f(mmo)
430 2691472896 : f(-m) = f(-mmo)*cp - sp*f(mmo)
431 : end do
432 :
433 271482624 : end subroutine ff
434 :
435 134281728 : subroutine pnm( ct, p )
436 : !----------------------------------------------------------------------------
437 : ! Purpose: normalized associated Legendre polynomial P_n^m
438 : ! Ref.: Richmond J.Atm.Ter.Phys. 1974
439 : ! Method:
440 : ! P_m^m = sqrt(1+1/2m)*si*P_m-1^m-1 m>0
441 : ! P_n^m = [cos*P_n-1^m - R_n-1^m*P_n-2^m ]/R_n^m n>m>=0
442 : ! dP/d phi = n*cos*P_n^m/sin-(2*n+1)*R_n^m*P_n-1^m/sin n>=m>=0
443 : ! R_n^m = sqrt[ (n^2-m^2)/(4n^2-1) ]
444 : !
445 : ! Author: A. Maute Nov 2003 am 11/18/03
446 : !----------------------------------------------------------------------------
447 :
448 : implicit none
449 :
450 : !-----------------------------------------------------------------------
451 : ! dummy arguments
452 : !-----------------------------------------------------------------------
453 : real(r8), intent(inout) :: ct ! cos(colat)
454 : real(r8), intent(inout) :: p(0:nm,0:mm)
455 :
456 : !-----------------------------------------------------------------------
457 : ! local variables
458 : !-----------------------------------------------------------------------
459 : integer :: mp, m, n
460 : real(r8) :: pm2, st
461 :
462 : ! ct = min( ct,.99_r8 ) ! cos(colat)
463 134281728 : st = sqrt( 1._r8 - ct*ct ) ! sin(colat)
464 :
465 134281728 : p(0,0) = 1._r8
466 2685634560 : do mp = 1,mmp ! m+1=1,mm+1
467 2551352832 : m = mp - 1
468 2551352832 : if( m >= 1 ) then
469 2417071104 : p(m,m) = pmopmmo(m)*p(m-1,m-1)*st
470 : end if
471 : pm2 = 0._r8
472 28199162880 : do n = mp,nm ! n=m+1,N
473 25513528320 : p(n,m) = (ct*p(n-1,m) - r(n-1,m)*pm2)/r(n,m)
474 28064881152 : pm2 = p(n-1,m)
475 : end do
476 : end do
477 :
478 134281728 : end subroutine pnm
479 :
480 1536 : subroutine prep_pnm
481 : !----------------------------------------------------------------------------
482 : ! Purpose: constant factors for normalized associated Legendre polynomial P_n^m
483 : ! Ref.: Richmond J.Atm.Ter.Phys. 1974
484 : !
485 : ! Method:
486 : ! PmoPmmo(m) = sqrt(1+1/2m)
487 : ! R_n^m = sqrt[ (n^2-m^2)/(4n^2-1) ]
488 : !
489 : ! Author: A. Maute Nov 2003 am 11/18/03
490 : !----------------------------------------------------------------------------
491 :
492 : implicit none
493 :
494 : !-----------------------------------------------------------------------
495 : ! local variables
496 : !-----------------------------------------------------------------------
497 : integer :: mp, m, n
498 : real(r8) :: xms, xns, den
499 :
500 30720 : do mp = 1, mmp ! m+1 = 1,mm+1
501 29184 : m = mp - 1
502 29184 : xms = m*m
503 29184 : if( mp /= 1 ) then
504 27648 : pmopmmo(m) = sqrt( 1._r8 + .5_r8/M )
505 : end if
506 351744 : do n = m,nm ! n = m,N
507 321024 : xns = n*n
508 321024 : den = max(4._r8*xns - 1._r8,1._r8)
509 350208 : r(n,m) = sqrt( (xns - xms)/den )
510 : end do
511 : end do
512 :
513 1536 : end subroutine prep_pnm
514 :
515 1536 : subroutine index_quiet
516 : !----------------------------------------------------------------------------
517 : ! Purpose: set up index for factors f_m(mlt),f_l(UT),f_-k(d) to
518 : ! describe the electric potential Phi for the empirical model
519 : !
520 : ! Method:
521 : ! Phi = sum_k sum_l sum_m sum_n [ A_klmn * P_n^m *f_m(mlt)*f_l(UT)*f_-k(d)]
522 : ! - since the electric potential is symmetric about the equator
523 : ! n+m odd terms are set zero resp. not used
524 : ! - in the summation for calculation Phi the index have the following
525 : ! range n=1,12 and m=-n,n, k=0,2 l=-2,2
526 : !
527 : ! Author: A. Maute Nov 2003 am 11/18/03
528 : !----------------------------------------------------------------------------
529 :
530 : implicit none
531 :
532 : !----------------------------------------------------------------------------
533 : ! ... local variables
534 : !----------------------------------------------------------------------------
535 : integer :: i, j, k, l, n, m
536 :
537 1536 : i = 0 ! initialize
538 1536 : j = 1
539 6144 : do k = 2,0,-1
540 29184 : do l = -2,2
541 23040 : if( k == 2 .and. abs(l) == 2 ) then
542 : cycle
543 : end if
544 264192 : do n = 1,12
545 9128448 : do m = -18,18
546 9105408 : if( abs(m) <= n ) then ! |m| < n
547 3354624 : if( (((n-m)/2)*2) == (n-m) ) then ! only n+m even
548 1797120 : if( n-abs(m) <= 9 ) then ! n-|m| <= 9 why?
549 1677312 : kf(i) = 2-k
550 1677312 : lf(i) = l
551 1677312 : nf(i) = n
552 1677312 : mf(i) = m
553 1677312 : jf(i) = j
554 1677312 : i = i + 1 ! counter
555 : end if
556 : end if
557 : end if
558 : end do ! m
559 : end do ! n
560 : end do ! l
561 : end do ! k
562 :
563 1536 : imax = i - 1
564 1536 : if(imax /= ni ) then ! check if imax == ni
565 0 : write(iulog,'(a19,i5,a18,i5)') 'index_quiet: imax= ',imax,' not equal to ni =',ni
566 0 : call endrun('index_quiet ERROR')
567 : end if
568 : if(debug.and.masterproc) write(iulog,*) 'imax=',imax
569 :
570 1536 : end subroutine index_quiet
571 :
572 1536 : subroutine read_acoef (efield_lflux_file, efield_hflux_file)
573 : !----------------------------------------------------------------------------
574 : ! Purpose:
575 : ! 1. read in coefficients A_klmn^lf for solar cycle minimum and
576 : ! A_klmn^hf for maximum
577 : ! 2. adjust S_a (f107d) such that if S_a<80 or S_a > 220 it has reasonable numbers
578 : ! S_aM = [atan{(S_a-65)^2/90^2}-a90]/[a180-a90]
579 : ! a90 = atan [(90-65)/90]^2
580 : ! a180 = atan [(180-65)/90]^2
581 : ! 3. inter/extrapolation of the coefficient to the actual flux which is
582 : ! given by the user
583 : ! A_klmn = S_aM [A_klmn^hf-A_klmn^lf]/90. + 2*A_klmn^lf-A_klmn^hf
584 : !
585 : ! Method:
586 : !
587 : ! Author: A. Maute Nov 2003 am 11/19/03
588 : !----------------------------------------------------------------------------
589 :
590 : use ioFileMod, only : getfil
591 : use units, only : getunit, freeunit
592 : use spmd_utils, only : mpicom, masterprocid, mpicom, mpi_real8
593 :
594 : character(len=*), intent(in) :: efield_lflux_file
595 : character(len=*), intent(in) :: efield_hflux_file
596 :
597 : integer :: ios,unit, ierr
598 : character (len=256):: locfn
599 : character(len=80) :: errmsg
600 :
601 1536 : if (masterproc) then
602 : !----------------------------------------------------------------------------
603 : ! get coefficients file for solar minimum:
604 : !----------------------------------------------------------------------------
605 2 : unit = getunit()
606 2 : call getfil( efield_lflux_file, locfn, 0 )
607 :
608 : !----------------------------------------------------------------------------
609 : ! open datafile with coefficients A_klnm
610 : !----------------------------------------------------------------------------
611 : if (debug) write(iulog,*) 'read_acoef: open file ',trim(locfn),' unit ',unit
612 : open(unit=unit,file=trim(locfn), &
613 2 : status = 'old',iostat = ios)
614 2 : if(ios.gt.0) then
615 0 : write(errmsg,*) 'read_acoef: error in opening coeff_lf file',' unit ',unit
616 0 : call endrun(errmsg)
617 : end if
618 :
619 : !----------------------------------------------------------------------------
620 : ! read datafile with coefficients A_klnm
621 : !----------------------------------------------------------------------------
622 : if (debug) write(iulog,*) 'read_acoef: read file ',trim(locfn),' unit ',unit
623 2 : read(unit,*,iostat = ios) a_lf
624 2 : if(ios.gt.0) then
625 0 : write(errmsg,*) 'read_acoef: error in reading coeff_lf file',' unit ',unit
626 0 : call endrun(errmsg)
627 : end if
628 :
629 : !----------------------------------------------------------------------------
630 : ! close & free unit
631 : !----------------------------------------------------------------------------
632 2 : close(unit)
633 2 : call freeunit(unit)
634 : if (debug) write(iulog,*) 'read_acoef: free unit ',unit
635 :
636 : !----------------------------------------------------------------------------
637 : ! get coefficients file for solar maximum:
638 : !----------------------------------------------------------------------------
639 2 : unit = getunit()
640 2 : call getfil( efield_hflux_file, locfn, 0 )
641 :
642 : !----------------------------------------------------------------------------
643 : ! open datafile with coefficients A_klnm
644 : !----------------------------------------------------------------------------
645 : if (debug) write(iulog,*) 'read_acoef: open file ',trim(locfn),' unit ',unit
646 : open(unit=unit,file=trim(locfn), &
647 2 : status = 'old',iostat = ios)
648 2 : if(ios.gt.0) then
649 0 : write(errmsg,*) 'read_acoef: error in opening coeff_hf file',' unit ',unit
650 0 : call endrun(errmsg)
651 : end if
652 :
653 : !----------------------------------------------------------------------------
654 : ! read datafile with coefficients A_klnm
655 : !----------------------------------------------------------------------------
656 : if (debug) write(iulog,*) 'read_acoef: read file ',trim(locfn)
657 2 : read(unit,*,iostat = ios) a_hf
658 2 : if(ios.gt.0) then
659 0 : write(errmsg,*) 'read_acoef: error in reading coeff_hf file',' unit ',unit
660 0 : call endrun(errmsg)
661 : end if
662 :
663 : !----------------------------------------------------------------------------
664 : ! close & free unit
665 : !----------------------------------------------------------------------------
666 2 : close(unit)
667 2 : call freeunit(unit)
668 : if (debug) write(iulog,*) 'read_acoef: free unit ',unit
669 : endif
670 :
671 1536 : call mpi_bcast (a_lf, ni+1, mpi_real8, masterprocid, mpicom, ierr)
672 1536 : call mpi_bcast (a_hf, ni+1, mpi_real8, masterprocid, mpicom, ierr)
673 :
674 1536 : end subroutine read_acoef
675 :
676 16128 : subroutine adj_S_a
677 : !----------------------------------------------------------------------------
678 : ! adjust S_a -> S_aM eqn.8-11 Scherliess draft
679 : !----------------------------------------------------------------------------
680 :
681 : implicit none
682 :
683 : !----------------------------------------------------------------------------
684 : ! local variables
685 : !----------------------------------------------------------------------------
686 : integer :: i
687 : real(r8) :: x2, y2, a90, a180, S_aM
688 :
689 16128 : x2 = 90._r8*90._r8
690 : y2 = (90._r8 - 65._r8)
691 : y2 = y2*y2
692 16128 : a90 = atan2(y2,x2)
693 : y2 = (180._r8 - 65._r8)
694 16128 : y2 = y2*y2
695 16128 : a180 = atan2(y2,x2)
696 : ! y2 = (S_a-65._r8)
697 16128 : y2 = (f107d - 65._r8)
698 16128 : y2 = y2*y2
699 16128 : S_aM = (atan2(y2,x2) - a90)/(a180 - a90)
700 16128 : S_aM = 90._r8*(1._r8 + S_aM)
701 : if(debug.and.masterproc) write(iulog,*) 'f107d=',f107d,' S_aM =',S_aM
702 :
703 : !----------------------------------------------------------------------------
704 : ! inter/extrapolate to S_a (f107d)
705 : !----------------------------------------------------------------------------
706 17627904 : do i = 0,ni ! eqn.8 Scherliess draft
707 17627904 : a_klnm(i) = S_aM*(a_hf(i) - a_lf(i))/90._r8 + 2._r8*a_lf(i) - a_hf(i)
708 : ! for testing like in original code
709 : ! a_klnm(i)=S_a*(a_hf(i)-a_lf(i))/90.+2.*a_lf(i)-a_hf(i)
710 : ! a_klnm(i)=f107d*(a_hf(i)-a_lf(i))/90.+2.*a_lf(i)-a_hf(i)
711 : end do
712 :
713 16128 : end subroutine adj_S_a
714 :
715 1536 : subroutine constants
716 : !----------------------------------------------------------------------------
717 : ! Purpose: set up constant values (e.g. magnetic grid, convertion
718 : ! constants etc)
719 : !
720 : ! Method:
721 : !
722 : ! Author: A. Maute Nov 2003 am 11/19/03
723 : !----------------------------------------------------------------------------
724 :
725 : !----------------------------------------------------------------------------
726 : ! local variables
727 : !----------------------------------------------------------------------------
728 : integer :: i,j
729 : real(r8) :: fac,lat
730 :
731 1536 : rtd = 180._r8/pi ! radians -> deg
732 1536 : dtr = pi/180._r8 ! deg -> radians
733 1536 : sqr2 = sqrt(2.e0_r8)
734 1536 : dy2rd = 2._r8*pi/dy2yr ! 2*pi/365.24 average year
735 :
736 : !----------------------------------------------------------------------------
737 : ! Set grid deltas:
738 : !----------------------------------------------------------------------------
739 1536 : dlatm = 180._r8/nmlat
740 1536 : dlonm = 360._r8/nmlon
741 :
742 : !----------------------------------------------------------------------------
743 : ! Set magnetic latitude array
744 : !----------------------------------------------------------------------------
745 141312 : do j = 0,nmlat
746 139776 : ylatm(j) = j*dlatm
747 139776 : lat = (ylatm(j) - 90._r8)*dtr
748 139776 : fac = cos(lat) ! sinIm = 2*sin(lam_m)/sqrt[4-3*cos^2(lam_m)]
749 139776 : fac = 4._r8 - 3._r8*fac*fac
750 139776 : fac = 2._r8/sqrt( fac )
751 141312 : sinIm_mag(j) = fac*sin( lat )
752 : end do
753 :
754 : !----------------------------------------------------------------------------
755 : ! Set magnetic longitude array
756 : !----------------------------------------------------------------------------
757 279552 : do i = 0,nmlon
758 279552 : ylonm(i) = i*dlonm
759 : end do
760 :
761 : !----------------------------------------------------------------------------
762 : ! find boundary index for high-lat potential
763 : !----------------------------------------------------------------------------
764 35328 : do j = 0,nmlat
765 35328 : nmlat_hgh = j
766 35328 : if( bnd_hgh <= ylatm(j) ) then
767 : exit
768 : end if
769 : end do
770 :
771 : !----------------------------------------------------------------------------
772 : ! find latitudinal shift
773 : !----------------------------------------------------------------------------
774 43008 : do j = 0,nmlat
775 43008 : ilat_sft = j
776 43008 : if( lat_sft <= ylatm(j) ) then
777 : exit
778 : end if
779 : end do
780 :
781 : !----------------------------------------------------------------------------
782 : ! find index for linear interpolation of ed2 at mag.equator
783 : ! use 12 deg - same as in TIEGCM
784 : !----------------------------------------------------------------------------
785 81408 : do j = 0,nmlat
786 81408 : lat = ylatm(j) - 90._r8
787 81408 : if( lat <= -12._r8 ) then
788 61440 : jmin = j
789 19968 : else if( lat > 12._r8 ) then
790 1536 : jmax = j
791 1536 : exit
792 : end if
793 : end do
794 :
795 1536 : end subroutine constants
796 :
797 1536 : subroutine prep_fk
798 : !----------------------------------------------------------------------------
799 : ! Purpose: set up constants factors for f_-k(day) used for empirical model
800 : ! to calculate the electric potential
801 : !
802 : ! Method:
803 : !
804 : ! Author: A. Maute Nov 2003 am 11/19/03
805 : !----------------------------------------------------------------------------
806 :
807 1536 : ft(1,0) = .75_r8*sqrt( 6.e0_r8 )/pi
808 1536 : ft(1,1) = 2.e0_r8*ft(1,0)
809 1536 : ft(1,2) = 1.e0_r8
810 1536 : ft(2,0) = ft(1,0)
811 1536 : ft(2,1) = -ft(1,1)
812 1536 : ft(2,2) = 1.e0_r8
813 1536 : ft(3,0) = ft(2,1)
814 1536 : ft(3,1) = 0._r8
815 1536 : ft(3,2) = 1.e0_r8
816 :
817 1536 : end subroutine prep_fk
818 :
819 134281728 : subroutine set_fkflfs( fk, fl, fs )
820 : !----------------------------------------------------------------------------
821 : ! Purpose: set f_-k(day) depending on seasonal flag used for empirical model
822 : ! to calculate the electric potential
823 : !
824 : ! Method:
825 : !
826 : ! Author: A. Maute Nov 2003 am 11/20/03
827 : !----------------------------------------------------------------------------
828 :
829 : !----------------------------------------------------------------------------
830 : ! ... dummy arguments
831 : !----------------------------------------------------------------------------
832 : real(r8), intent(out) :: &
833 : fk(0:2), & ! f_-k(day)
834 : fl(-2:2), & ! f_l(ut)
835 : fs(2) ! f_s(f10.7)
836 : !----------------------------------------------------------------------------
837 : ! local variables
838 : !----------------------------------------------------------------------------
839 : integer :: lp
840 : real(r8) :: ang
841 : real(r8) :: lon_ut
842 :
843 : !----------------------------------------------------------------------------
844 : ! f_-k(day)
845 : ! use factors for iseasav == 0 - Scherliess had iseasav as an input parameter
846 : !----------------------------------------------------------------------------
847 134281728 : lp = iseasav
848 : if( iseasav == 0 ) then
849 134281728 : ang = (day + 9._r8)*dy2rd
850 134281728 : fk(0) = sqr2*cos( 2._r8*ang )
851 134281728 : fk(1) = sqr2*cos( ang )
852 134281728 : fk(2) = 1._r8
853 : else if( iseasav >= 1 .and. iseasav <= 3 ) then
854 : fk(0) = ft(lp,0)
855 : fk(1) = ft(lp,1)
856 : fk(2) = ft(lp,2)
857 : else if( iseasav == 4 ) then
858 : fk(0) =0._r8
859 : fk(1) =0._r8
860 : fk(2) =1._r8
861 : end if
862 :
863 : !----------------------------------------------------------------------------
864 : ! f_l(ut)
865 : !----------------------------------------------------------------------------
866 134281728 : lon_ut = 15._r8*ut ! 15.*mlt - xmlon + 69.
867 134281728 : call ff( lon_ut, 2, fl )
868 : if( iutav ) then ! UT-averaging
869 :
870 : ang = fl(0)
871 : fl(:) = 0._r8
872 : fl(0) = ang
873 :
874 : end if
875 :
876 : !----------------------------------------------------------------------------
877 : ! f_s(f10.7) only fs(1) used
878 : !----------------------------------------------------------------------------
879 134281728 : fs(1) = 1._r8
880 : ! fs(2) = S_a
881 134281728 : fs(2) = f107d
882 :
883 134281728 : end subroutine set_fkflfs
884 :
885 134281728 : subroutine efield_mid( mlat, mlon, pot )
886 : !----------------------------------------------------------------------------
887 : ! Purpose: calculate the electric potential for low and
888 : ! midlatitudes from an empirical model (Scherliess 1999)
889 : !
890 : ! Method:
891 : !
892 : ! Author: A. Maute Nov 2003 am 11/20/03
893 : !----------------------------------------------------------------------------
894 :
895 : !----------------------------------------------------------------------------
896 : ! ... dummy arguments
897 : !----------------------------------------------------------------------------
898 : real(r8), intent(in) :: mlat, mlon
899 : real(r8), intent(out) :: pot ! electric potential (V)
900 :
901 : !----------------------------------------------------------------------------
902 : ! local variables
903 : !----------------------------------------------------------------------------
904 : integer :: i, mp, np, nn
905 : real(r8) :: mod_mlat, ct, x
906 : real(r8) :: fk(0:2) ! f_-k(day)
907 : real(r8) :: fl(-2:2) ! f_l(ut)
908 : real(r8) :: fs(2) ! f_s(f10.7)
909 : real(r8) :: f(-18:18)
910 : real(r8) :: p(0:nm,0:mm) ! P_n^m spherical harmonics
911 :
912 134281728 : pot = 0._r8 ! initialize
913 :
914 134281728 : mod_mlat = mlat
915 134281728 : if( abs(mlat) <= 0.5_r8 ) then
916 2919168 : mod_mlat = 0.5_r8 ! avoid geomag.equator
917 : end if
918 :
919 : !----------------------------------------------------------------------------
920 : ! set f_-k, f_l, f_s depending on seasonal flag
921 : !----------------------------------------------------------------------------
922 134281728 : call set_fkflfs( fk, fl, fs )
923 :
924 : !----------------------------------------------------------------------------
925 : ! spherical harmonics
926 : !----------------------------------------------------------------------------
927 134281728 : ct = cos( (90._r8 - mod_mlat)*dtr ) ! magnetic colatitude
928 134281728 : call pnm( ct, p ) ! calculate P_n^m
929 134281728 : call ff( mlon, 18, f ) ! calculate f_m (phi) why 18 if N=12
930 :
931 >14676*10^7 : do i = 0,imax
932 >14663*10^7 : mp = mf(i)
933 >14663*10^7 : np = nf(i)
934 >14663*10^7 : nn = abs(mp) ! P_n^m = P_n^-m
935 >14663*10^7 : x = a_klnm(i)* fl(lf(i)) * fk(kf(i)) * fs(jf(i))
936 >14676*10^7 : pot = pot + x*f(mp)*p(np,nn)
937 : end do
938 :
939 134281728 : end subroutine efield_mid
940 :
941 16128 : subroutine pot_latsmo( pot, idlat ) ! pots == pot_highlats
942 : !----------------------------------------------------------------------------
943 : ! Purpose: smoothing in latitude of potential
944 : !
945 : ! Method: weighted smoothing in latitude
946 : ! assume regular grid spacing
947 : !
948 : ! Author: A. Maute Nov 2003 am 11/20/03
949 : !----------------------------------------------------------------------------
950 :
951 : !----------------------------------------------------------------------------
952 : ! ... dummy arguments
953 : !----------------------------------------------------------------------------
954 : integer, intent(in) :: idlat
955 : real(r8), intent(inout) :: pot(0:nmlon,0:nmlat)
956 :
957 : !----------------------------------------------------------------------------
958 : ! local variables
959 : !----------------------------------------------------------------------------
960 : integer :: ilat, id
961 : real(r8) :: wgt, del
962 32256 : real(r8) :: w(-idlat:idlat)
963 : ! real(r8) :: pot_smo(0:nmlat) ! temp array for smooth. potential
964 32256 : real(r8) :: pot_smo(0:nmlon,0:nmlat_hgh) ! temp array for smooth. potential
965 :
966 : !----------------------------------------------------------------------------
967 : ! weighting factors (regular grid spacing)
968 : !----------------------------------------------------------------------------
969 16128 : wgt = 0._r8
970 96768 : do id = -idlat,idlat
971 80640 : del = abs(id)*dlatm ! delta lat_m
972 80640 : w(id) = 1._r8/(del + 1._r8)
973 96768 : wgt = wgt + w(id)
974 : end do
975 16128 : wgt = 1._r8/wgt
976 :
977 : !$omp parallel do private(ilat)
978 322560 : do ilat = idlat,nmlat_hgh-idlat
979 55786752 : pot_smo(:,ilat) = matmul( pot(:,ilat-idlat:ilat+idlat),w )*wgt
980 : end do
981 :
982 322560 : do ilat = idlat,nmlat_hgh-idlat
983 56077056 : pot(:,ilat) = pot_smo(:,ilat)
984 55786752 : pot(:,nmlat-ilat) = pot_smo(:,ilat)
985 : end do
986 :
987 16128 : end subroutine pot_latsmo
988 :
989 16128 : subroutine pot_latsmo2( pot, idlat )
990 : !----------------------------------------------------------------------------
991 : ! Purpose: smoothing in latitude of potential
992 : !
993 : ! Method: weighted smoothing in latitude
994 : ! assume regular grid spacing
995 : !
996 : ! Author: A. Maute Nov 2003 am 11/20/03
997 : !----------------------------------------------------------------------------
998 :
999 : !----------------------------------------------------------------------------
1000 : ! ... dummy arguments
1001 : !----------------------------------------------------------------------------
1002 : integer, intent(in) :: idlat
1003 : real(r8), intent(inout) :: pot(0:nmlon,0:nmlat)
1004 :
1005 : !----------------------------------------------------------------------------
1006 : ! local variables
1007 : !----------------------------------------------------------------------------
1008 : integer :: ilat, id
1009 : real(r8) :: wgt, del
1010 32256 : real(r8) :: w(-idlat:idlat)
1011 : ! real(r8) :: pot_smo(0:nmlat) ! temp array for smooth. potential
1012 : real(r8) :: pot_smo(0:nmlon,0:nmlath) ! temp array for smooth. potential
1013 :
1014 : !----------------------------------------------------------------------------
1015 : ! weighting factors (regular grid spacing)
1016 : !----------------------------------------------------------------------------
1017 16128 : wgt = 0._r8
1018 96768 : do id = -idlat,idlat
1019 80640 : del = abs(id)*dlatm ! delta lat_m
1020 80640 : w(id) = 1._r8/(del + 1._r8)
1021 96768 : wgt = wgt + w(id)
1022 : end do
1023 16128 : wgt = 1._r8/wgt
1024 :
1025 : ! do ilon = 0,nmlon
1026 : ! do ilat = idlat,nmlath-idlat ! ilat = 5:175
1027 : ! pot_smo(ilat) = 0._r8
1028 : ! do id = -idlat,idlat ! org. was degree now grid points
1029 : ! pot_smo(ilat) = pot_smo(ilat) + w(id)*pot(ilon,ilat+id)
1030 : ! end do
1031 : ! pot_smo(ilat) = pot_smo(ilat)*wgt
1032 : ! end do
1033 : ! pot(ilon,idlat:nmlath-idlat) = pot_smo(idlat:nmlath-idlat) ! copy back into pot
1034 : ! end do
1035 :
1036 : !$omp parallel do private(ilat)
1037 693504 : do ilat = idlat,nmlath-idlat
1038 123298560 : pot_smo(:,ilat) = matmul( pot(:,ilat-idlat:ilat+idlat),w )*wgt
1039 : end do
1040 :
1041 693504 : do ilat = idlat,nmlath-idlat
1042 123975936 : pot(:,ilat) = pot_smo(:,ilat)
1043 : end do
1044 :
1045 16128 : end subroutine pot_latsmo2
1046 :
1047 16128 : subroutine pot_lonsmo( pot, idlon )
1048 : !----------------------------------------------------------------------------
1049 : ! Purpose: smoothing in longitude of potential
1050 : !
1051 : ! Method: weighted smoothing in longitude
1052 : ! assume regular grid spacing
1053 : !
1054 : ! Author: A. Maute Nov 2003 am 11/20/03
1055 : !----------------------------------------------------------------------------
1056 :
1057 : !----------------------------------------------------------------------------
1058 : ! ... dummy arguments
1059 : !----------------------------------------------------------------------------
1060 : integer, intent(in) :: idlon
1061 : real(r8), intent(inout) :: pot(0:nmlon,0:nmlat)
1062 :
1063 : !----------------------------------------------------------------------------
1064 : ! local variables
1065 : !----------------------------------------------------------------------------
1066 : integer :: ilon, ilat, id
1067 : real(r8) :: wgt, del
1068 32256 : real(r8) :: w(-idlon:idlon)
1069 32256 : real(r8) :: tmp(-idlon:nmlon+idlon) ! temp array for smooth. potential
1070 :
1071 : !----------------------------------------------------------------------------
1072 : ! weighting factors (regular grid spacing)
1073 : !----------------------------------------------------------------------------
1074 16128 : wgt = 0._r8
1075 129024 : do id = -idlon,idlon
1076 112896 : del = abs(id)*dlonm ! delta lon_m
1077 112896 : w(id) = 1._r8/(del + 1._r8)
1078 129024 : wgt = wgt + w(id)
1079 : end do
1080 16128 : wgt = 1._r8/wgt
1081 :
1082 : !----------------------------------------------------------------------------
1083 : ! averaging
1084 : !----------------------------------------------------------------------------
1085 : ! do ilon = 0,nmlon
1086 : ! do ilat = 0,nmlath
1087 : ! pot_smo(ilat) = 0._r8
1088 : ! do id = -idlon,idlon ! org. was degree now grid points
1089 : ! iabs = ilon + id
1090 : ! if( iabs > nmlon ) then
1091 : ! iabs = iabs - nmlon ! test if wrap around
1092 : ! end if
1093 : ! if( iabs < 0 ) then
1094 : ! iabs = iabs + nmlon ! test if wrap around
1095 : ! end if
1096 : ! pot_smo(ilat) = pot_smo(ilat) + w(id)*pot(iabs,ilat)
1097 : ! end do
1098 : ! pot_smo(ilat) = pot_smo(ilat)*wgt
1099 : ! pot(ilon,ilat) = pot_smo(ilat) ! copy back into pot
1100 : ! pot(ilon,nmlat-ilat) = pot_smo(ilat) ! copy back into pot
1101 : ! end do
1102 : ! end do
1103 :
1104 : !$omp parallel do private(ilat,ilon,tmp)
1105 758016 : do ilat = 0,nmlath
1106 135023616 : tmp(0:nmlon) = pot(0:nmlon,ilat)
1107 2967552 : tmp(-idlon:-1) = pot(nmlon-idlon:nmlon-1,ilat)
1108 3709440 : tmp(nmlon+1:nmlon+idlon) = pot(1:idlon,ilat)
1109 135039744 : do ilon = 0,nmlon
1110 1074253824 : pot(ilon,ilat) = dot_product( tmp(ilon-idlon:ilon+idlon),w )*wgt
1111 135023616 : pot(ilon,nmlat-ilat) = pot(ilon,ilat)
1112 : end do
1113 : end do
1114 :
1115 16128 : end subroutine pot_lonsmo
1116 :
1117 16128 : subroutine highlat_getbnd( ihlat_bnd )
1118 : !----------------------------------------------------------------------------
1119 : ! Purpose: calculate the height latitude bounday index ihl_bnd
1120 : !
1121 : ! Method:
1122 : ! 1. calculate E field from high-latitude model
1123 : ! boundary is set where the total electric field exceeds
1124 : ! 0.015 V/m (corresp. approx. to 300 m/s)
1125 : ! 2. moved halfways to 54 deg not necessarily equatorwards as in the
1126 : ! original comment from L. Scherliess- or?
1127 : !
1128 : ! Author: A. Maute Nov 2003 am 11/20/03
1129 : !----------------------------------------------------------------------------
1130 :
1131 : !----------------------------------------------------------------------------
1132 : ! ... dummy arguments
1133 : !----------------------------------------------------------------------------
1134 : integer, intent(out) :: ihlat_bnd(0:nmlon)
1135 :
1136 : !----------------------------------------------------------------------------
1137 : ! local variables
1138 : !----------------------------------------------------------------------------
1139 : integer :: ilon, ilat, ilat_sft_rvs
1140 : real(r8) :: e_tot
1141 : real(r8) :: e_max, ef_max
1142 :
1143 : ! first find absolute max E-field magnitude
1144 16128 : e_max = 0._r8
1145 2935296 : do ilon = 0,nmlon
1146 72995328 : do ilat = nmlat_hgh+1,0,-1 ! lat. loop moving torwards pole
1147 70060032 : e_tot = sqrt( ed1(ilon,ilat)**2 + ed2(ilon,ilat)**2 )
1148 72979200 : if (e_tot > e_max) then
1149 677376 : e_max = e_tot
1150 : end if
1151 : end do
1152 : end do
1153 :
1154 : ! Set E-field strength used to find the transition boundary:
1155 : ! when Kp is low the max E-field strength from Heelis is less than 0.015 V/m
1156 : ! for such cases set ef_max to half the max E-field strength
1157 16128 : ef_max = min(0.015_r8,e_max*0.5_r8)
1158 :
1159 16128 : ilat_sft_rvs = nmlath - ilat_sft ! pole =0, equ=90
1160 2935296 : do ilon = 0,nmlon
1161 46480896 : do ilat = nmlat_hgh+1,0,-1 ! lat. loop moving torwards pole
1162 46464768 : e_tot = sqrt( ed1(ilon,ilat)**2 + ed2(ilon,ilat)**2 )
1163 46464768 : if( abs(e_tot) >= ef_max ) then ! e-filed > limit -> boundary
1164 2919168 : ihlat_bnd(ilon) = ilat - (ilat - ilat_sft_rvs)*.5_r8 ! shift boundary to lat_sft (54deg)
1165 2919168 : exit
1166 : end if
1167 : end do
1168 : end do
1169 :
1170 16128 : end subroutine highlat_getbnd
1171 :
1172 16128 : subroutine bnd_sinus( ihlat_bnd, itrans_width )
1173 : !----------------------------------------------------------------------------
1174 : ! Purpose:
1175 : ! 1. adjust high latitude boundary (ihlat_bnd) sinusoidally
1176 : ! 2. width of transition zone from midlatitude potential to high latitude
1177 : ! potential (itrans_width)
1178 : !
1179 : ! Method:
1180 : ! 1.adjust boundary sinusoidally
1181 : ! max. wave number to be represented nmax_sin
1182 : ! RHS(mi) = Sum_phi Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi)*hlat_bnd(phi)
1183 : ! U(mi,mk) = Sum_phi Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi) *
1184 : ! Sum_(mk=-nmax_sin)^_(mk=nmax_sin) f_mk(phi)
1185 : ! single values decomposition of U
1186 : ! solving U*LSG = RHS
1187 : ! calculating hlat_bnd:
1188 : ! hlat_bnd = Sum_(mi=-nmax_sin)^_(mi=nmax_sin) f_mi(phi)*LSG(mi)
1189 : !
1190 : ! 2. width of transition zone from midlatitude potential to high latitude
1191 : ! potential
1192 : ! trans_width(phi)=8.-2.*cos(phi)
1193 : !
1194 : ! Author: A. Maute Nov 2003 am 11/20/03
1195 : !----------------------------------------------------------------------------
1196 :
1197 : external DGESV ! LAPACK routine to solve matrix eq
1198 :
1199 : !----------------------------------------------------------------------------
1200 : ! ... dummy arguments
1201 : !----------------------------------------------------------------------------
1202 : integer, intent(inout) :: ihlat_bnd(0:nmlon) ! loaction of boundary
1203 : integer, intent(out) :: itrans_width(0:nmlon) ! width of transition zone
1204 :
1205 : !----------------------------------------------------------------------------
1206 : ! local variables
1207 : !----------------------------------------------------------------------------
1208 : integer, parameter :: nmax_a = 2*nmax_sin+1 ! absolute array length
1209 : integer, parameter :: ishf = nmax_sin+1
1210 : integer :: ilon, i, i1, j, bnd
1211 : real(r8) :: sum
1212 : real(r8) :: rhs(nmax_a)
1213 : real(r8) :: lsg(nmax_a)
1214 : real(r8) :: u(nmax_a,nmax_a)
1215 : real(r8) :: v(nmax_a,nmax_a)
1216 : real(r8) :: w(nmax_a,nmax_a)
1217 : real(r8) :: f(-nmax_sin:nmax_sin,0:nmlon)
1218 :
1219 : real(r8) :: x(nmax_a)
1220 : integer :: ipiv(nmax_a), info
1221 :
1222 : character(len=120) :: msg
1223 :
1224 : !----------------------------------------------------------------------------
1225 : ! Sinusoidal Boundary calculation
1226 : !----------------------------------------------------------------------------
1227 16128 : rhs(:) = 0._r8
1228 16128 : lsg(:) = 0._r8
1229 16128 : u(:,:) = 0._r8
1230 16128 : v(:,:) = 0._r8
1231 16128 : w(:,:) = 0._r8
1232 16128 : ipiv(:) = 0
1233 :
1234 2935296 : do ilon = 0,nmlon ! long.
1235 2919168 : bnd = nmlath - ihlat_bnd(ilon) ! switch from pole=0 to pole =90
1236 2919168 : call ff( ylonm(ilon), nmax_sin, f(-nmax_sin,ilon) )
1237 17531136 : do i = -nmax_sin,nmax_sin
1238 14595840 : i1 = i + ishf
1239 14595840 : rhs(i1) = rhs(i1) + f(i,ilon) * bnd
1240 : ! if (debug) write(iulog,*) 'rhs ',ilon,i1,bnd,f(i,ilon),rhs(i+ishf)
1241 90494208 : do j = -nmax_sin,nmax_sin
1242 87575040 : u(i1,j+ishf) = u(i1,j+ishf) + f(i,ilon)*f(j,ilon)
1243 : ! if (debug) write(iulog,*) 'u ',ilon,i1,j+ishf,u(i+ishf,j+ishf)
1244 : end do
1245 : end do
1246 : end do
1247 : !
1248 16128 : x(:) = rhs(:)
1249 16128 : call DGESV( nmax_a, 1, u, nmax_a, ipiv, x, nmax_a, info)
1250 16128 : if (info/=0) then
1251 0 : write(msg,'(a,i4)') 'bnd_sinus -- LAPACK DGESV return error code: ',info
1252 0 : if (masterproc) write(iulog,*) trim(msg)
1253 0 : call endrun(trim(msg))
1254 : end if
1255 16128 : lsg(:) = x(:)
1256 : !
1257 2935296 : do ilon = 0,nmlon ! long.
1258 17515008 : sum = dot_product( lsg(-nmax_sin+ishf:nmax_sin+ishf),f(-nmax_sin:nmax_sin,ilon) )
1259 2919168 : ihlat_bnd(ilon) = nmlath - int( sum + .5_r8 ) ! closest point
1260 2935296 : itrans_width(ilon) = int( 8._r8 - 2._r8*cos( ylonm(ilon)*dtr ) + .5_r8 )/dlatm ! 6 to 10 deg.
1261 : end do
1262 : ! if (debug) write(iulog,"('bnd_sinus: ihlat_bnd=',/,(12i6))") ihlat_bnd
1263 : ! if (debug) write(iulog,"('bnd_sinus: itrans_width=',/,(12i6))") itrans_width
1264 :
1265 16128 : end subroutine bnd_sinus
1266 :
1267 16128 : subroutine highlat_adjust( pot_highlats, pot_highlat, pot_midlat, ihlat_bnd )
1268 : !----------------------------------------------------------------------------
1269 : ! Purpose: Adjust mid/low latitude electric potential and high latitude
1270 : ! potential such that there are continous across the mid to high
1271 : ! latitude boundary
1272 : !
1273 : ! Method:
1274 : ! 1. integrate Phi_low/mid(phi,bnd) along the boundary mid to high latitude
1275 : ! 2. integrate Phi_high(phi,bnd) along the boundary mid to high latitude
1276 : ! 3. adjust Phi_high by delta =
1277 : ! Int_phi Phi_high(phi,bnd) d phi/360. - Int_phi Phi_low/mid(phi,bnd) d phi/360.
1278 : !
1279 : ! Author: A. Maute Nov 2003 am 11/21/03
1280 : !----------------------------------------------------------------------------
1281 :
1282 : !----------------------------------------------------------------------------
1283 : ! ... dummy arguments
1284 : !----------------------------------------------------------------------------
1285 : integer, intent(in) :: ihlat_bnd(0:nmlon) ! boundary mid to high latitude
1286 : real(r8), intent(in) :: pot_midlat(0:nmlon,0:nmlat) ! low/mid latitude potentail
1287 : real(r8), intent(inout) :: pot_highlat(0:nmlon,0:nmlat) ! high_lat potential
1288 : real(r8), intent(inout) :: pot_highlats(0:nmlon,0:nmlat) ! high_lat potential! smoothed high_lat potential
1289 :
1290 : !----------------------------------------------------------------------------
1291 : ! local:
1292 : !----------------------------------------------------------------------------
1293 : integer :: ilon, ilat, ilatS, ibnd60, ibnd_hl
1294 : real(r8) :: pot60, pot_hl, del
1295 :
1296 : !----------------------------------------------------------------------------
1297 : ! 1. integrate Phi_low/mid(phi,bnd) along the boundary mid to high latitude
1298 : ! 2. integrate Phi_high(phi,bnd) along the boundary mid to high latitude
1299 : !----------------------------------------------------------------------------
1300 16128 : pot60 = 0._r8
1301 16128 : pot_hl = 0._r8
1302 2919168 : do ilon = 1,nmlon ! long. ! bnd -> eq to pole 0:90
1303 2903040 : ibnd60 = nmlat - ihlat_bnd(ilon) ! 0:180 pole to pole
1304 2903040 : ibnd_hl = ihlat_bnd(ilon) ! colatitude
1305 2903040 : pot60 = pot60 + pot_midlat(ilon,ibnd60)
1306 2919168 : pot_hl = pot_hl + pot_highlats(ilon,ibnd_hl)
1307 : end do
1308 16128 : pot60 = pot60/(nmlon)
1309 16128 : pot_hl = pot_hl/(nmlon)
1310 :
1311 : if (debug.and.masterproc) write(iulog,*) 'Mid-Latitude Boundary Potential =',pot60
1312 : if (debug.and.masterproc) write(iulog,*) 'High-Latitude Boundary Potential=',pot_hl
1313 :
1314 : !----------------------------------------------------------------------------
1315 : ! 3. adjust Phi_high by delta =
1316 : ! Int_phi Phi_high(phi,bnd) d phi/360. - Int_phi Phi_low/mid(phi,bnd) d phi/360.
1317 : !----------------------------------------------------------------------------
1318 16128 : del = pot_hl - pot60
1319 :
1320 : !$omp parallel do private(ilat,ilon,ilats)
1321 387072 : do ilat = 0,nmlat_hgh ! colatitude
1322 370944 : ilats = nmlat - ilat
1323 67527936 : do ilon = 0,nmlon
1324 67140864 : pot_highlat(ilon,ilat) = pot_highlat(ilon,ilat) - del
1325 67140864 : pot_highlat(ilon,ilats) = pot_highlat(ilon,ilats) - del
1326 67140864 : pot_highlats(ilon,ilat) = pot_highlats(ilon,ilat) - del
1327 67511808 : pot_highlats(ilon,ilats) = pot_highlats(ilon,ilats) - del
1328 : end do
1329 : end do
1330 :
1331 16128 : end subroutine highlat_adjust
1332 :
1333 16128 : subroutine interp_poten( pot_highlats, pot_highlat, pot_midlat, &
1334 : ihlat_bnd, itrans_width )
1335 : !----------------------------------------------------------------------------
1336 : ! Purpose: construct a smooth global electric potential field
1337 : !
1338 : ! Method: construct one global potential field
1339 : ! 1. low/mid latitude: |lam| < bnd-trans_width
1340 : ! Phi(phi,lam) = Phi_low(phi,lam)
1341 : ! 2. high latitude: |lam| > bnd+trans_width
1342 : ! Phi(phi,lam) = Phi_hl(phi,lam)
1343 : ! 3. transition zone: bnd-trans_width <= lam <= bnd+trans_width
1344 : ! a. interpolate between high and low/midlatitude potential
1345 : ! Phi*(phi,lam) = 1/15*[ 5/(2*trans_width) * {Phi_low(phi,bnd-trans_width)*
1346 : ! [-lam+bnd+trans_width] + Phi_hl(phi,bnd+trans_width)*
1347 : ! [lam-bnd+trans_width]} + 10/(2*trans_width) {Phi_low(phi,lam)*
1348 : ! [-lam+bnd+trans_width] + Phi_hl(phi,lam)*
1349 : ! [lam-bnd+trans_width]}]
1350 : ! b. Interpolate between just calculated Potential and the high latitude
1351 : ! potential in a 3 degree zone poleward of the boundary:
1352 : ! bnd+trans_width < lam <= bnd+trans_width+ 3 deg
1353 : ! Phi(phi,lam) = 1/3 { [3-(lam-bnd-trans_width)]* Phi*(phi,lam) +
1354 : ! [lam-bnd-trans_width)]* Phi_hl*(phi,lam) }
1355 : !
1356 : ! Author: A. Maute Nov 2003 am 11/21/03
1357 : !----------------------------------------------------------------------------
1358 :
1359 : !----------------------------------------------------------------------------
1360 : ! ... dummy arguments
1361 : !----------------------------------------------------------------------------
1362 : integer, intent(in) :: ihlat_bnd(0:nmlon)
1363 : integer, intent(in) :: itrans_width(0:nmlon)
1364 : real(r8), intent(in) :: pot_highlats(0:nmlon,0:nmlat)
1365 : real(r8), intent(in) :: pot_highlat(0:nmlon,0:nmlat)
1366 : real(r8), intent(in) :: pot_midlat(0:nmlon,0:nmlat)
1367 :
1368 : !----------------------------------------------------------------------------
1369 : ! local variables
1370 : !----------------------------------------------------------------------------
1371 : real(r8), parameter :: fac = 1._r8/3._r8
1372 : integer :: ilon, ilat
1373 : integer :: ibnd, tw, hb1, hb2, lat_ind
1374 : integer :: j1, j2
1375 : real(r8) :: a, b, b1, b2
1376 : real(r8) :: wrk1, wrk2
1377 :
1378 : !$omp parallel do private(ilat,ilon,ibnd,tw)
1379 2935296 : do ilon = 0,nmlon
1380 2919168 : ibnd = ihlat_bnd(ilon) ! high latitude boundary index
1381 2919168 : tw = itrans_width(ilon) ! width of transition zone (index)
1382 : !----------------------------------------------------------------------------
1383 : ! 1. low/mid latitude: |lam| < bnd-trans_width
1384 : ! Phi(phi,lam) = Phi_low(phi,lam)
1385 : !----------------------------------------------------------------------------
1386 86155776 : do ilat = 0,nmlath-(ibnd+tw+1)
1387 83236608 : potent(ilon,nmlath+ilat) = pot_midlat(ilon,nmlath+ilat)
1388 86155776 : potent(ilon,nmlath-ilat) = pot_midlat(ilon,nmlath+ilat)
1389 : end do
1390 : !----------------------------------------------------------------------------
1391 : ! 2. high latitude: |lam| > bnd+trans_width
1392 : ! Phi(phi,lam) = Phi_hl(phi,lam)
1393 : !----------------------------------------------------------------------------
1394 28836864 : do ilat = 0,ibnd-tw-1
1395 25901568 : potent(ilon,ilat) = pot_highlats(ilon,nmlat-ilat)
1396 28820736 : potent(ilon,nmlat-ilat) = pot_highlats(ilon,nmlat-ilat)
1397 : end do
1398 : end do
1399 : !----------------------------------------------------------------------------
1400 : ! 3. transition zone: bnd-trans_width <= lam <= bnd+trans_width
1401 : !----------------------------------------------------------------------------
1402 : ! a. interpolate between high and low/midlatitude potential
1403 : ! update only southern hemisphere (northern hemisphere is copied
1404 : ! after smoothing)
1405 : !----------------------------------------------------------------------------
1406 : !$omp parallel do private(ilat,ilon,ibnd,tw,a,b,b1,b2,hb1,hb2,lat_ind,j1,j2,wrk1,wrk2)
1407 2935296 : do ilon = 0,nmlon
1408 2919168 : ibnd = ihlat_bnd(ilon) ! high latitude boundary index
1409 2919168 : tw = itrans_width(ilon) ! width of transition zone (index)
1410 2919168 : a = 1._r8/(2._r8*tw)
1411 2919168 : b1 = (nmlath - ibnd + tw)*a
1412 2919168 : b2 = (nmlath - ibnd - tw)*a
1413 2919168 : hb1 = nmlath - (ibnd + tw)
1414 2919168 : j1 = nmlath - hb1
1415 2919168 : hb2 = nmlath - (ibnd - tw)
1416 2919168 : j2 = nmlath - hb2
1417 2919168 : wrk1 = pot_midlat(ilon,j1)
1418 2919168 : wrk2 = pot_highlats(ilon,j2)
1419 : ! if(debug.and.masterproc) write(iulog,*) 'pot_all ',ilon,hb1,hb2,nmlath -ibnd,tw
1420 28062720 : do ilat = ibnd-tw,ibnd+tw
1421 25143552 : lat_ind = nmlath - ilat
1422 25143552 : potent(ilon,ilat) = &
1423 : fac*((wrk1 + 2._r8*pot_midlat(ilon,ilat))*(b1 - a*lat_ind) &
1424 25143552 : + (wrk2 + 2._r8*pot_highlats(ilon,ilat))*(a*lat_ind - b2))
1425 28062720 : potent(ilon,nmlat-ilat) = potent(ilon,ilat)
1426 : end do
1427 : !----------------------------------------------------------------------------
1428 : ! b. Interpolate between just calculated Potential and the high latitude
1429 : ! potential in a 3 degree zone poleward of the boundary
1430 : !----------------------------------------------------------------------------
1431 28836864 : do ilat = hb2+1,nmlath
1432 25901568 : a = max( 3._r8/dlatm - (ilat - hb2 - 1),0._r8 )
1433 25901568 : b = 3._r8/dlatm - a
1434 25901568 : potent(ilon,nmlath-ilat) = (a*potent(ilon,nmlath-ilat) &
1435 25901568 : + b*pot_highlat(ilon,nmlath-ilat))/3._r8*dlatm
1436 28820736 : potent(ilon,nmlath+ilat) = potent(ilon,nmlath-ilat)
1437 : end do
1438 : end do
1439 :
1440 16128 : end subroutine interp_poten
1441 :
1442 32256 : subroutine DerivPotential
1443 : !-----------------------------------------------------------------------
1444 : ! Purpose: calulates the electric field [V/m] from the electric potential
1445 : !
1446 : ! Method: Richmond [1995] eqn 5.9-5.10
1447 : ! ed1(:,:) = Ed1 = - 1/[R cos lam_m] d PHI/d phi_m
1448 : ! ed2(:,:) = Ed2 = 1/R d PHI/d lam_m /sin I_m
1449 : ! R = R_e + h_r we assume a reference height of 130 km which is also
1450 : ! used in the TIEGCM code
1451 : !
1452 : ! Author: A. Maute Dec 2003 am 12/16/03
1453 : !-----------------------------------------------------------------------
1454 :
1455 : integer :: i, j
1456 : real(r8) :: coslm, r, fac, wrk
1457 : real(r8) :: wrk1d(0:nmlon)
1458 :
1459 32256 : r = r_e + h_r ! earth radius + reference height [m]
1460 : !-----------------------------------------------------------------------
1461 : ! ed2= Ed2 is the equatorward/downward component of the electric field, at all
1462 : ! geomagnetic grid points (central differencing)
1463 : !-----------------------------------------------------------------------
1464 32256 : fac = .5_r8/(dlatm*dtr*r)
1465 : !$omp parallel do private(j, i, wrk )
1466 1451520 : do j = 1,nmlath-1 ! southern hemisphere
1467 1419264 : wrk = fac/sinIm_mag(j)
1468 258338304 : do i = 0,nmlon
1469 258306048 : ed2(i,j) = (potent(i,j+1) - potent(i,j-1))*wrk
1470 : end do
1471 : end do
1472 :
1473 : !$omp parallel do private(j, i, wrk )
1474 1451520 : do j = nmlath+1,nmlat-1 ! northern hemisphere
1475 1419264 : wrk = fac/sinIm_mag(j)
1476 258338304 : do i = 0,nmlon
1477 258306048 : ed2(i,j) = (potent(i,j+1) - potent(i,j-1))*wrk
1478 : end do
1479 : end do
1480 :
1481 : !-----------------------------------------------------------------------
1482 : ! Interpolate of ed2 between between -12 <= lam_m <= 12 degrees:
1483 : !-----------------------------------------------------------------------
1484 5870592 : wrk1d(:) = ed2(:,jmax) - ed2(:,jmin)
1485 : !$omp parallel do private(j, i, fac)
1486 419328 : do j = jmin+1,jmax-1
1487 387072 : fac = (ylatm(j) - ylatm(jmin))/(ylatm(jmax) - ylatm(jmin))
1488 70479360 : do i = 0,nmlon
1489 70447104 : ed2(i,j) = ed2(i,jmin) + fac*wrk1d(i)
1490 : end do
1491 : end do
1492 :
1493 : !-----------------------------------------------------------------------
1494 : ! ed1= Ed1 is the zonal component of the electric field, at all
1495 : ! geomagnetic grid points (central differencing)
1496 : !-----------------------------------------------------------------------
1497 32256 : fac = .5_r8/(dlonm*dtr*r)
1498 : !$omp parallel do private(j, i, wrk, coslm )
1499 2903040 : do j = 1,nmlat-1
1500 2870784 : coslm = ylatm(j) - 90._r8
1501 2870784 : coslm = cos( coslm*dtr )
1502 2870784 : wrk = fac/coslm
1503 516741120 : do i = 1,nmlon-1
1504 516741120 : ed1(i,j) = -(potent(i+1,j) - potent(i-1,j))*wrk
1505 : end do
1506 2870784 : i = 0
1507 2870784 : ed1(i,j) = -(potent(i+1,j) - potent(nmlon-1,j))*wrk
1508 2903040 : ed1(nmlon,j) = ed1(i,j)
1509 : end do
1510 :
1511 : !-----------------------------------------------------------------------
1512 : ! Poles:
1513 : !-----------------------------------------------------------------------
1514 : !$omp parallel do private(i)
1515 5870592 : do i = 0,nmlon
1516 5838336 : ed1(i,0) = .25_r8*(ed1(i,1) - ed1(ip2f(i),1) + ed2(ip1f(i),1) - ed2(ip3f(i),1))
1517 : ed1(i,nmlat) = .25_r8*(ed1(i,nmlat-1) - ed1(ip2f(i),nmlat-1) &
1518 5838336 : + ed2(ip1f(i),nmlat-1) - ed2(ip3f(i),nmlat-1))
1519 5838336 : ed2(i,0) = .25_r8*(ed2(i,1) - ed2(ip2f(i),1) - ed1(ip1f(i),1) + ed1(ip3f(i),1))
1520 : ed2(i,nmlat) = .25_r8*(ed2(i,nmlat-1) - ed2(ip2f(i),nmlat-1) &
1521 5870592 : - ed1(ip1f(i),nmlat-1) + ed1(ip3f(i),nmlat-1))
1522 : end do
1523 :
1524 32256 : end subroutine DerivPotential
1525 :
1526 : end module efield
|
