Line data Source code
1 :
2 : module exbdrift
3 : !----------------------------------------------------------------------
4 : ! description: calculates ExB drift velocities UI,VI,WI
5 : ! uses the electric field which is calculated in module efield
6 : ! on a regular magnetic grid (MLT deg/ mag. latitude)
7 : !
8 : ! 0. initilize called before time-loop (exbdrift_init)
9 : ! every timestep and every processor
10 : ! 1. map from magn. grid to geographic grid WACCM (map_mag2geo)
11 : ! 2. rotate e-field (rot_efield)
12 : ! 3. calculate ExB drift velocities ui,vi,wi [m/s] (iondrift)
13 : !
14 : ! input ed1,ed2, & ! zonal/meridional elect. field [V/m]
15 : ! nmlon,nmlat, & ! dimension of mag. grid
16 : ! dlatm,dlonm, & ! grid spacing of mag. grid
17 : ! ylatm,ylonm ! magnetic MLT deg./latitudes (deg)
18 : !
19 : ! ExB electromagnetic drift velocity [m/s] (east,north,upward)
20 : ! These are output to the physics buffer.
21 : ! real(r8) :: ui,vi,wi
22 : !
23 : ! notes:
24 : ! - assume regular magnetic grid for the e-field (0:360 MLT/0:90) -> mapping
25 : !
26 : ! Author: A. Maute Dec 2003
27 : ! B. Foster adding physics buffer Feb, 2004.
28 : !----------------------------------------------------------------------
29 :
30 : use shr_kind_mod, only: r8 => shr_kind_r8
31 : use physconst, only: pi
32 :
33 : use ppgrid , only: pcols, pver
34 : use cam_logfile, only: iulog
35 :
36 : use efield, only : & ! inputs from efield module
37 : ed1, ed2, & ! global zonal/meridional elect. field [V/m]
38 : potent, & ! electric potential [V]
39 : nmlon, nmlat, & ! dimension of mag. grid
40 : dlatm, dlonm, & ! grid spacing of mag. grid
41 : ylatm, ylonm ! magnetic longitudes,latitudes (deg) (0:nmlat),(0:nmlon)
42 : use apex, only : apex_subsol, apex_magloctm
43 : use cam_history, only: outfld, addfld, add_default, horiz_only ! for history saves
44 : use shr_assert_mod, only: shr_assert_in_domain
45 : use phys_control, only: phys_getopts
46 :
47 : implicit none
48 :
49 : private
50 :
51 : save
52 :
53 : !----------------------------------------------------------------------
54 : ! Public interfaces:
55 : !----------------------------------------------------------------------
56 : public :: exbdrift_init
57 : public :: exbdrift_register ! register drift velocities with pbuf
58 : public :: exbdrift_ion_vels ! updates empirical ion drift velocities (in pbuf)
59 :
60 : !----------------------------------------------------------------------
61 : ! Indices to drift velocities in physics buffer:
62 : !----------------------------------------------------------------------
63 : integer :: ndx_ui, ndx_vi, ndx_wi
64 : real(r8), parameter :: rtd = 180._r8/pi ! radians to degrees
65 : real(r8), parameter :: hr2d = 360._r8/24._r8
66 :
67 : logical :: state_debug_checks = .false.
68 :
69 : contains
70 :
71 74496 : subroutine exbdrift_init( empirical_ion_vels )
72 : !-----------------------------------------------------------------------
73 : ! Purpose: Prepare fields for histories.
74 : !
75 : ! Method:
76 : !
77 : ! Author: A. Maute Dec 2003 am 12/30/03
78 : !-----------------------------------------------------------------------
79 :
80 : use phys_control, only: phys_getopts
81 :
82 : logical, intent(in) :: empirical_ion_vels
83 :
84 : logical :: history_waccm
85 :
86 :
87 1536 : call phys_getopts(history_waccm_out=history_waccm, state_debug_checks_out=state_debug_checks)
88 :
89 : !-----------------------------------------------------------------------
90 : ! Add mag field output to master field list:
91 : !-----------------------------------------------------------------------
92 1536 : call addfld('EF_EAST', horiz_only,'I','V/m', 'eastward electric field')
93 1536 : call addfld('EF_NORTH', horiz_only,'I','V/m', 'northward electric field')
94 1536 : call addfld('EF_UP', horiz_only,'I','V/m', 'upward electric field')
95 1536 : call addfld('EF1_MAP', horiz_only,'I','V/m', 'map. mag. eastward ef')
96 1536 : call addfld('EF2_MAP', horiz_only,'I','V/m', 'map. mag. northward ef')
97 1536 : call addfld('EPOTEN', horiz_only,'I','V', 'Electric Potential')
98 : !-----------------------------------------------------------------------
99 : ! Write these fields to WACCM history by default:
100 : !-----------------------------------------------------------------------
101 1536 : if (history_waccm .and. empirical_ion_vels) then
102 1536 : call add_default ('EPOTEN ' , 1, ' ')
103 : end if
104 :
105 1536 : end subroutine exbdrift_init
106 :
107 1536 : subroutine exbdrift_register
108 :
109 : use physics_buffer, only : pbuf_add_field, dtype_r8
110 :
111 : !-----------------------------------------------------------------------
112 : ! Register drift velocity outputs with physics buffer:
113 : !
114 : ! Ion velocities are 2d fields in pbuf (no vertical dimension),
115 : ! so fdim = mdim = ldim = 1.
116 : ! Indices are saved as module data above.
117 : !-----------------------------------------------------------------------
118 :
119 1536 : call pbuf_add_field("UI", "global", dtype_r8, (/pcols,pver/), ndx_ui)
120 1536 : call pbuf_add_field("VI", "global", dtype_r8, (/pcols,pver/), ndx_vi)
121 1536 : call pbuf_add_field("WI", "global", dtype_r8, (/pcols,pver/), ndx_wi)
122 :
123 1536 : end subroutine exbdrift_register
124 :
125 72960 : subroutine exbdrift_ion_vels( lchnk, ncol, pbuf)
126 1536 : use physics_buffer, only : physics_buffer_desc
127 : !-----------------------------------------------------------------------
128 : ! Purpose: calculate ion drift velocities [m/s]
129 : !
130 : ! Method: v = E x B/B^2
131 : !
132 : ! Author: A. Maute Dec 2003 am 12/30/03
133 : !-----------------------------------------------------------------------
134 :
135 : !-----------------------------------------------------------------------
136 : ! ... dummy arguments
137 : !-----------------------------------------------------------------------
138 : integer, intent(in) :: ncol ! np. of atmospheric columns
139 : integer, intent(in) :: lchnk ! chunk identifier
140 :
141 : type(physics_buffer_desc), pointer :: pbuf(:)
142 :
143 : !-----------------------------------------------------------------------
144 : ! ... local variables
145 : !-----------------------------------------------------------------------
146 145920 : real(r8) :: ed1_geo(ncol), & ! electric field on geographic grid [V/m]
147 145920 : ed2_geo(ncol), &
148 145920 : epot_geo(ncol) ! electric potential on geographic grid
149 145920 : real(r8) :: elfld(3,ncol) ! electric field in geog. direction on geographic grid [V/m]
150 145920 : real(r8) :: mlt(ncol) ! mag.local time of WACCM geo. grid point
151 :
152 : !-----------------------------------------------------------------------
153 : ! calculate the magn. local time of WACCM geogr. grid points
154 : !-----------------------------------------------------------------------
155 72960 : call cal_mlt( mlt, lchnk, ncol )
156 : !-----------------------------------------------------------------------
157 : ! map the electric field from regular mag.grid to WACCM grid
158 : !-----------------------------------------------------------------------
159 72960 : call map_mag2geo( mlt, lchnk, ncol, ed1_geo, ed2_geo, epot_geo )
160 : !-----------------------------------------------------------------------
161 : ! rotate the electric field in geographic direction
162 : !-----------------------------------------------------------------------
163 72960 : call rot_efield( lchnk, ncol, ed1_geo, ed2_geo, elfld )
164 : !-----------------------------------------------------------------------
165 : ! calculate ExB drift velocities ui,vi,wi [m/s]
166 : !-----------------------------------------------------------------------
167 72960 : call iondrift( lchnk, ncol, elfld, pbuf)
168 :
169 72960 : end subroutine exbdrift_ion_vels
170 :
171 72960 : subroutine map_mag2geo( mlt, lchnk, pcol, ed1_geo, ed2_geo, epot_geo )
172 : !-----------------------------------------------------------------------
173 : ! Purpose: map electric field from regular magnetic grid
174 : ! to WACCM physics geographic grid
175 : !
176 : ! Method: bilinear interpolation
177 : ! assumptions: magnetic grid regular (MLT[deg]=0:360 deg/lat=0:180 deg)
178 : ! ed1_geo(:,:) = Ed1_geo <- mapping of - 1/[R cos lam_m] d PHI/d phi_m
179 : ! ed2_geo(:,:) = Ed2_geo <- mapping of 1/R d PHI/d lam_m/ sinIm
180 : !
181 : ! Author: A.Maute Dec 2003 am 12/17/03
182 : !-----------------------------------------------------------------------
183 :
184 72960 : use mo_apex, only: alatm ! apex mag latitude at each geographic grid point (radians)
185 :
186 : !-----------------------------------------------------------------------
187 : ! dummy arguments
188 : !-----------------------------------------------------------------------
189 : integer, intent(in) :: pcol ! np. of atmospheric columns
190 : integer, intent(in) :: lchnk ! chunk identifier
191 : real(r8), intent(in) :: mlt(pcol) ! mag.local time of WACCM geo. grid point
192 : real(r8), intent(out) :: &
193 : ed1_geo(pcol), & ! electric field on geog. grid
194 : ed2_geo(pcol), & ! electric field on geog. grid
195 : epot_geo(pcol) ! electric potential on geog grid
196 :
197 : !-----------------------------------------------------------------------
198 : ! local variables
199 : !-----------------------------------------------------------------------
200 : integer :: i, iphi1, iphi2, ilam1, ilam2
201 : real(r8) :: t, u, collat, collon
202 :
203 72960 : if (state_debug_checks) then
204 72960 : call shr_assert_in_domain(ed1, is_nan=.false., varname="ed1", msg="NaN found in exbdrift::map_mag2geo")
205 72960 : call shr_assert_in_domain(ed2, is_nan=.false., varname="ed2", msg="NaN found in exbdrift::map_mag2geo")
206 72960 : call shr_assert_in_domain(potent, is_nan=.false., varname="potent", msg="NaN found in exbdrift::map_mag2geo")
207 : end if
208 :
209 1123584 : do i = 1,pcol
210 1050624 : collat = alatm(i,lchnk)*rtd ! mag lats (deg) [-90:90]
211 1050624 : collon = mlt(i)*hr2d ! mlt (deg) [0:360]
212 1050624 : iphi1 = int(collon/dlonm) ! indices of the 4 surrounding points (regular mag. grid
213 1050624 : iphi2 = iphi1 + 1 ! with dlonm/dlatm deg grid spacing)
214 1050624 : ilam1 = int( (collat + 90._r8)/dlatm )
215 1050624 : ilam2 = ilam1 + 1
216 1050624 : if(iphi1 == nmlon ) then ! boundaries
217 0 : iphi1 = nmlon -1
218 0 : iphi2 = nmlon
219 : end if
220 1050624 : if(ilam1 == nmlat ) then
221 0 : ilam1 = nmlat-1
222 0 : ilam2 = nmlat
223 : end if
224 1050624 : if(collon == 360._r8) then
225 0 : collon= 0._r8
226 0 : iphi1 = 0
227 0 : iphi2 = 1
228 : end if
229 :
230 1050624 : t = (collon - ylonm(iphi1))/(ylonm(iphi2) - ylonm(iphi1))
231 1050624 : u = (collat + 90._r8 - ylatm(ilam1))/(ylatm(ilam2) - ylatm(ilam1))
232 : ed1_geo(i) = (1._r8 - t)*(1._r8 - u)*ed1(iphi1,ilam1) + &
233 : t*(1._r8 - u)*ed1(iphi2,ilam1) + &
234 : t*u* ed1(iphi2,ilam2) + &
235 1050624 : (1._r8 - t)*u*ed1(iphi1,ilam2)
236 : ed2_geo(i) = (1._r8 - t)*(1._r8 - u)*ed2(iphi1,ilam1) + &
237 : t*(1._r8 - u)*ed2(iphi2,ilam1) + &
238 : t*u* ed2(iphi2,ilam2) + &
239 1050624 : (1._r8 - t)*u*ed2(iphi1,ilam2)
240 : epot_geo(i)= (1._r8 - t)*(1._r8 - u)*potent(iphi1,ilam1) + &
241 : t*(1._r8 - u)*potent(iphi2,ilam1) + &
242 : t*u* potent(iphi2,ilam2) + &
243 1123584 : (1._r8 - t)*u*potent(iphi1,ilam2)
244 :
245 : end do ! i = 1,pcol
246 :
247 72960 : if (state_debug_checks) then
248 72960 : call shr_assert_in_domain(ed1_geo, is_nan=.false., varname="ed1_geo", msg="NaN found in exbdrift::map_mag2geo")
249 72960 : call shr_assert_in_domain(ed2_geo, is_nan=.false., varname="ed2_geo", msg="NaN found in exbdrift::map_mag2geo")
250 72960 : call shr_assert_in_domain(epot_geo, is_nan=.false., varname="epot__geo", msg="NaN found in exbdrift::map_mag2geo")
251 : endif
252 :
253 72960 : call outfld( 'EF1_MAP', ed1_geo, pcol, lchnk)
254 72960 : call outfld( 'EF2_MAP', ed2_geo, pcol, lchnk)
255 72960 : call outfld( 'EPOTEN', epot_geo, pcol, lchnk)
256 :
257 72960 : end subroutine map_mag2geo
258 :
259 72960 : subroutine rot_efield( lchnk, pcol, ed1_geo, ed2_geo, elfld )
260 : !-----------------------------------------------------------------------
261 : ! Purpose: rotate the electric field to get geographic east and westward direction
262 : !
263 : ! Method: Richmond: J Geomag. Geoelectr. [1995] eqn. (4.5)
264 : ! rotation of the electric field to get geog. eastward and downward/equatorward
265 : ! E(k) = d_1(k)*ed1_geo + d_2(k)*ed2_geo for k = 1,3
266 : !
267 : ! Author: A. Maute Dec 2003 am 12/17/03
268 : !-----------------------------------------------------------------------
269 :
270 : use mo_apex, only: d1vec, d2vec ! base vectors (3,pcols,begchunk:endchunk)
271 :
272 : !-----------------------------------------------------------------------
273 : ! dummy arguments
274 : !-----------------------------------------------------------------------
275 : integer, intent(in) :: pcol ! np. of atmospheric columns
276 : integer, intent(in) :: lchnk ! chunk identifier
277 : real(r8), intent(in ) :: ed1_geo(pcol),ed2_geo(pcol) ! electric field on geog. grid
278 : real(r8), intent(out) :: elfld(3,pcol) ! electric field on geog. grid geog. direction
279 :
280 : !-----------------------------------------------------------------------
281 : ! local
282 : !-----------------------------------------------------------------------
283 : integer :: i, k
284 :
285 1123584 : do i = 1,pcol
286 4275456 : do k = 1,3
287 6303744 : elfld(k,i) = ed1_geo(i)*d1vec(k,i,lchnk) + &
288 10506240 : ed2_geo(i)*d2vec(k,i,lchnk)
289 : end do
290 : end do
291 :
292 1123584 : call outfld( 'EF_EAST', elfld(1,:), pcol, lchnk)
293 1123584 : call outfld( 'EF_NORTH', elfld(2,:), pcol, lchnk)
294 1123584 : call outfld( 'EF_UP', elfld(3,:), pcol, lchnk)
295 :
296 72960 : end subroutine rot_efield
297 :
298 145920 : subroutine iondrift( lchnk, pcol, elfld, pbuf)
299 : !-----------------------------------------------------------------------
300 : ! Purpose: calculate ion drift velocity
301 : !
302 : ! Method: v_i = ExB/B^2
303 : ! for high altitudes where collisionfrequency mue_in << omega_i and
304 : ! mue_en,vert << omega_e
305 : ! v_i,vert = v_e,vert = v_E
306 : ! B magnetic field component from apex code
307 : ! bnorth northward gauss
308 : ! beast eastward gauss
309 : ! bdown downward gauss -> for upward component -bdown
310 : ! bmag magnitude gauss
311 : !
312 : ! Author: A. Maute Dec 2003 am 12/17/03
313 : !-----------------------------------------------------------------------
314 :
315 : use mo_apex, only: beast, bnorth, bdown, bmag ! component of B-field, |B| [Gauss]
316 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field
317 :
318 :
319 : !-----------------------------------------------------------------------
320 : ! dummy arguments
321 : !-----------------------------------------------------------------------
322 : integer, intent(in) :: pcol ! np. of atmospheric columns
323 : integer, intent(in) :: lchnk ! chunk id
324 : real(r8), intent(in) :: elfld(3,pcol) ! electric field on geog. grid geog. direction
325 : !-----------------------------------------------------------------------
326 : ! Ion velocities are saved in physics buffer:
327 : !-----------------------------------------------------------------------
328 :
329 : type(physics_buffer_desc), pointer :: pbuf(:)
330 :
331 : !-----------------------------------------------------------------------
332 : ! local variables
333 : !-----------------------------------------------------------------------
334 : integer :: i
335 : real(r8) :: fac
336 72960 : real(r8), pointer :: ui(:,:), vi(:,:), wi(:,:)
337 :
338 : !-----------------------------------------------------------------------
339 : ! Set local pointers to respective fields in pbuf:
340 : !-----------------------------------------------------------------------
341 72960 : call pbuf_get_field(pbuf, ndx_ui, ui )
342 72960 : call pbuf_get_field(pbuf, ndx_vi, vi )
343 72960 : call pbuf_get_field(pbuf, ndx_wi, wi )
344 :
345 1123584 : do i = 1,pcol ! number of columns in each chunk
346 : !-----------------------------------------------------------------------
347 : ! Magnitude of magnetic field bmag [nT] is from apex module.
348 : !-----------------------------------------------------------------------
349 1050624 : fac = 1.e9_r8/bmag(i,lchnk)**2 ! nT to T (nT in denominator)
350 74594304 : ui(i,:) = -(elfld(2,i)*bdown(i,lchnk) + elfld(3,i)*bnorth(i,lchnk))
351 74594304 : vi(i,:) = elfld(3,i)*beast(i,lchnk) + elfld(1,i)*bdown(i,lchnk)
352 74594304 : wi(i,:) = elfld(1,i)*bnorth(i,lchnk) - elfld(2,i)*beast(i,lchnk)
353 74594304 : ui(i,:) = ui(i,:)*fac
354 74594304 : vi(i,:) = vi(i,:)*fac
355 74667264 : wi(i,:) = wi(i,:)*fac
356 : end do ! i = 1,pcol
357 :
358 72960 : if (state_debug_checks) then
359 72960 : call shr_assert_in_domain(ui(:pcol,:), is_nan=.false., varname="ui", msg="NaN found in exbdrift::iondrift")
360 72960 : call shr_assert_in_domain(vi(:pcol,:), is_nan=.false., varname="vi", msg="NaN found in exbdrift::iondrift")
361 72960 : call shr_assert_in_domain(wi(:pcol,:), is_nan=.false., varname="wi", msg="NaN found in exbdrift::iondrift")
362 : endif
363 :
364 : #ifdef SW_DEBUG
365 : if( lchnk == 25 ) then
366 : write(iulog,*) ' '
367 : write(iulog,*) '---------------------------------------'
368 : write(iulog,*) 'iondrift: elfld @ lchnk,i = ',lchnk,' 14'
369 : write(iulog,'(1p,3g15.7)') elfld(:,14)
370 : write(iulog,*) 'iondrift: bdown,bnorth,beast,bmag @ lchnk,i = ',lchnk,' 14'
371 : write(iulog,'(1p,4g15.7)') bdown(14,lchnk), bnorth(14,lchnk), beast(14,lchnk), bmag(14,lchnk)
372 : write(iulog,*) '---------------------------------------'
373 : write(iulog,*) ' '
374 : end if
375 : #endif
376 :
377 145920 : end subroutine iondrift
378 :
379 72960 : subroutine cal_mlt( mlt, lchnk, ncol )
380 :
381 : !-------------------------------------------------------------------------------
382 : ! Purpose: calculate the magnetic local time of WACCM geog. point
383 : !
384 : ! Method: using the location of the geomagnetic dipole north pole,
385 : ! the subsolar point location and the apex longitude of the
386 : ! geographic WACCM point the magn. local time can be calculated
387 : ! subroutines from Roy Barnes HAO Feb. 2004
388 : !
389 : ! Author: A. Maute Feb 2004
390 : !-------------------------------------------------------------------------------
391 :
392 : use mo_apex, only: &
393 : alonm, & ! apex mag longitude at each geographic grid point (radians)
394 72960 : colatp, & ! geocentric colatitude of geomagnetic dipole north pole (deg)
395 : elonp ! East longitude of geomagnetic dipole north pole (deg)
396 : use time_manager, only : get_curr_calday, get_curr_date
397 : use mo_apex, only : geomag_year
398 : !-------------------------------------------------------------------------------
399 : ! ... dummy arguments
400 : !-------------------------------------------------------------------------------
401 : integer, intent(in) :: ncol ! np. of atmospheric columns
402 : integer, intent(in) :: lchnk ! chunk identifier
403 : real(r8), intent(out) :: mlt(ncol) ! mag.local time of WACCM geo. grid point
404 :
405 : !-------------------------------------------------------------------------------
406 : ! local arguments
407 : !-------------------------------------------------------------------------------
408 : integer :: i
409 : integer :: iyear, iday, ihr, imn, imo, iday_m, tod ! time of day [s]
410 : real(r8) :: collonm, sbsllat, sbsllon, sec
411 :
412 : !-------------------------------------------------------------------------------
413 : ! get current calendar day of year & date components
414 : ! valid at end of current timestep
415 : !-------------------------------------------------------------------------------
416 72960 : call get_curr_date (iyear,imo,iday_m,tod) ! year, time of day [sec]
417 72960 : iyear = int(geomag_year)
418 72960 : iday = get_curr_calday() ! day of year
419 :
420 72960 : ihr = tod/3600
421 72960 : imn = mod( tod,3600 )/60
422 72960 : sec = tod - 60*(ihr*60 + imn)
423 :
424 : !-------------------------------------------------------------------------------
425 : ! find subsolar geographic latitude and longitude
426 : !-------------------------------------------------------------------------------
427 72960 : call apex_subsol( iyear, iday, ihr, imn, sec, sbsllat, sbsllon )
428 :
429 : !-------------------------------------------------------------------------------
430 : ! computes magnetic local time
431 : !-------------------------------------------------------------------------------
432 1123584 : do i = 1,ncol
433 1050624 : collonm = alonm(i,lchnk)*rtd ! mag lons (deg)
434 1123584 : call apex_magloctm( collonm, sbsllat, sbsllon, colatp, elonp, mlt(i) )
435 : end do
436 :
437 72960 : end subroutine cal_mlt
438 :
439 : end module exbdrift
|
