Line data Source code
1 :
2 : module iondrag
3 : !-------------------------------------------------------------------------------
4 : ! Purpose:
5 : ! Calculate ion drag tendency and apply to horizontal velocities.
6 : ! Also calculate joule heating tendency and apply to neutral temperature.
7 : !
8 : ! Subroutines:
9 : ! iondrag_init (initialize module)
10 : ! iondrag_calc (calculate ion drag tensors)
11 : ! iondrag_tend (ion drag tendency)
12 : ! qjoule_tend (joule heating)
13 : !
14 : ! Calling sequence:
15 : ! inti
16 : ! iondrag_init
17 : ! tphysac
18 : ! iondrag_calc
19 : ! iondrag_tend
20 : ! qjoule_tend
21 : !
22 : ! Dependencies:
23 : ! Magnetic field from apex module
24 : ! ExB ion drifts from exbdrift module
25 : !
26 : ! Author:
27 : ! B. Foster Feb, 2004
28 : !
29 : !-------------------------------------------------------------------------------
30 :
31 : use shr_kind_mod ,only: r8 => shr_kind_r8
32 : use ppgrid ,only: pcols, pver, begchunk, endchunk
33 : use cam_history ,only: addfld, add_default, outfld, horiz_only
34 : use physics_types,only: physics_state, physics_ptend, physics_ptend_init
35 :
36 : use physics_buffer, only : pbuf_get_index, physics_buffer_desc, pbuf_get_field
37 : use perf_mod ,only: t_startf, t_stopf
38 : use cam_logfile ,only: iulog
39 :
40 : use interpolate_data, only: lininterp
41 : use spmd_utils, only: masterproc
42 :
43 : use phys_control, only: waccmx_is
44 : use cam_abortutils,only: endrun
45 : use time_manager, only: is_first_step
46 :
47 : implicit none
48 :
49 : save
50 :
51 : private ! Make default type private to the module
52 :
53 : !-------------------------------------------------------------------------------
54 : ! Public interfaces:
55 : !-------------------------------------------------------------------------------
56 : public :: iondrag_register ! Register variables in pbuf physics buffer
57 : public :: iondrag_init ! Initialization
58 : public :: iondrag_calc ! ion drag tensors lxx,lyy,lxy,lyx
59 : public :: iondrag_readnl
60 : public :: iondrag_timestep_init
61 : public :: iondrag_inidat
62 : public :: do_waccm_ions
63 :
64 : interface iondrag_calc
65 : module procedure iondrag_calc_ions
66 : module procedure iondrag_calc_ghg
67 : end interface
68 :
69 : !-------------------------------------------------------------------------------
70 : ! Private data
71 : !-------------------------------------------------------------------------------
72 :
73 : ! Namelist variables
74 : character(len=256) :: efield_lflux_file = 'coeff_lflux.dat'
75 : character(len=256) :: efield_hflux_file = 'coeff_hflux.dat'
76 : real(r8) :: efield_potential_max = huge(1._r8) ! max cross cap potential kV
77 : logical :: empirical_ion_velocities = .true.
78 :
79 : real(r8),parameter :: amu = 1.6605387e-27_r8 ! atomic mass unit (kg)
80 :
81 : integer :: ntop_lev = 1
82 : integer :: nbot_lev = 0
83 : integer :: id_xo2, id_xo1 ! indices to tn and major sp
84 : integer :: id_o2p, id_op, id_nop ! indices to ions
85 : integer :: id_elec, id_n
86 :
87 : !Physics buffer indices
88 : integer :: PedConduct_idx = 0
89 : integer :: HallConduct_idx = 0
90 : integer :: ui_idx = 0 ! index to zonal drift from edynamo
91 : integer :: vi_idx = 0 ! index to meridional drift from edynamo
92 : integer :: wi_idx = 0 ! index to vertical drift from edynamo
93 : integer :: indxTe = -1 ! pbuf index for electron temperature
94 : integer :: indxTi = -1 ! pbuf index for ion temperature
95 :
96 : logical :: xo2_slvd, xo1_slvd, o2p_slvd, op_slvd, nop_slvd
97 :
98 : real(r8) :: rmass_op ! mass of O+ (g/mole)
99 : real(r8) :: rmass_o2p ! mass of O2+ (g/mole)
100 : real(r8) :: rmass_nop ! mass of NO+ (g/mole)
101 : real(r8) :: rmass_o1 ! mass of O (g/mole)
102 : real(r8) :: rmass_o2 ! mass of O2 (g/mole)
103 : real(r8) :: rmass_n2 ! mass of N2 (g/mole)
104 :
105 : !-------------------------------------------------------------------------------
106 : ! Inverted masses (for multiply in loops rather than divide):
107 : !-------------------------------------------------------------------------------
108 : real(r8) :: rmi_o1
109 : real(r8) :: rmi_o2
110 : real(r8) :: rmi_n2
111 : real(r8) :: rmi_op
112 : real(r8) :: rmi_o2p
113 : real(r8) :: rmi_nop
114 : real(r8) :: rmi_op_kg
115 : real(r8) :: rmi_o2p_kg
116 : real(r8) :: rmi_nop_kg
117 :
118 : ! GHG
119 : !-------------------------------------------------------------------------
120 :
121 : ! Private data
122 : integer, parameter :: plevtiod = 97
123 :
124 : real(r8) alamxx(plevtiod)
125 : real(r8) alamxy(plevtiod)
126 :
127 : real(r8) pshtiod(plevtiod) ! TIME pressure scale height
128 : real(r8) pshccm(pver) ! CCM pressure scale height
129 :
130 : real(r8) alamxxi(pver) ! alamxx interpolated to waccm grid
131 : real(r8) alamxyi(pver) ! alamxy interpoalted to waccm grid
132 :
133 : logical doiodrg
134 : logical, protected :: do_waccm_ions = .false.
135 :
136 : !
137 : ! Data statement for ALAMXX
138 : !
139 : data alamxx / &
140 : 0.13902E-17_r8, 0.22222E-17_r8, 0.34700E-17_r8, 0.53680E-17_r8, 0.83647E-17_r8, &
141 : 0.13035E-16_r8, 0.20254E-16_r8, 0.31415E-16_r8, 0.48944E-16_r8, 0.75871E-16_r8, &
142 : 0.11584E-15_r8, 0.17389E-15_r8, 0.25786E-15_r8, 0.37994E-15_r8, 0.58088E-15_r8, &
143 : 0.95179E-15_r8, 0.19052E-14_r8, 0.47869E-14_r8, 0.14284E-13_r8, 0.45584E-13_r8, &
144 : 0.14756E-12_r8, 0.48154E-12_r8, 0.14844E-11_r8, 0.39209E-11_r8, 0.83886E-11_r8, &
145 : 0.14213E-10_r8, 0.20304E-10_r8, 0.27449E-10_r8, 0.39276E-10_r8, 0.59044E-10_r8, &
146 : 0.83683E-10_r8, 0.11377E-09_r8, 0.14655E-09_r8, 0.19059E-09_r8, 0.28338E-09_r8, &
147 : 0.46326E-09_r8, 0.73966E-09_r8, 0.11785E-08_r8, 0.18789E-08_r8, 0.31037E-08_r8, &
148 : 0.53919E-08_r8, 0.97251E-08_r8, 0.17868E-07_r8, 0.33041E-07_r8, 0.61265E-07_r8, &
149 : 0.11406E-06_r8, 0.20912E-06_r8, 0.39426E-06_r8, 0.76691E-06_r8, 0.15113E-05_r8, &
150 : 0.29545E-05_r8, 0.55644E-05_r8, 0.97208E-05_r8, 0.16733E-04_r8, 0.28101E-04_r8, &
151 : 0.36946E-04_r8, 0.44277E-04_r8, 0.50982E-04_r8, 0.57526E-04_r8, 0.64190E-04_r8, &
152 : 0.71471E-04_r8, 0.80311E-04_r8, 0.96121E-04_r8, 0.11356E-03_r8, 0.14131E-03_r8, &
153 : 0.18695E-03_r8, 0.26058E-03_r8, 0.36900E-03_r8, 0.50812E-03_r8, 0.66171E-03_r8, &
154 : 0.80763E-03_r8, 0.92583E-03_r8, 0.10038E-02_r8, 0.10382E-02_r8, 0.10333E-02_r8, &
155 : 0.99732E-03_r8, 0.93994E-03_r8, 0.86984E-03_r8, 0.79384E-03_r8, 0.71691E-03_r8, &
156 : 0.64237E-03_r8, 0.57224E-03_r8, 0.50761E-03_r8, 0.44894E-03_r8, 0.39624E-03_r8, &
157 : 0.34929E-03_r8, 0.30767E-03_r8, 0.27089E-03_r8, 0.23845E-03_r8, 0.20985E-03_r8, &
158 : 0.18462E-03_r8, 0.16233E-03_r8, 0.14260E-03_r8, 0.12510E-03_r8, 0.10955E-03_r8, &
159 : 0.95699E-04_r8, 0.83347E-04_r8/
160 :
161 : !
162 : ! Data statement for ALAMXY
163 : !
164 : data alamxy / &
165 : 0.74471E-24_r8, 0.22662E-23_r8, 0.69004E-23_r8, 0.20345E-22_r8, 0.58465E-22_r8, &
166 : 0.16542E-21_r8, 0.46240E-21_r8, 0.12795E-20_r8, 0.35226E-20_r8, 0.96664E-20_r8, &
167 : 0.26650E-19_r8, 0.76791E-19_r8, 0.25710E-18_r8, 0.10897E-17_r8, 0.56593E-17_r8, &
168 : 0.30990E-16_r8, 0.16792E-15_r8, 0.85438E-15_r8, 0.40830E-14_r8, 0.18350E-13_r8, &
169 : 0.79062E-13_r8, 0.33578E-12_r8, 0.13348E-11_r8, 0.45311E-11_r8, 0.12443E-10_r8, &
170 : 0.27052E-10_r8, 0.49598E-10_r8, 0.86072E-10_r8, 0.15807E-09_r8, 0.30480E-09_r8, &
171 : 0.55333E-09_r8, 0.96125E-09_r8, 0.15757E-08_r8, 0.25896E-08_r8, 0.48209E-08_r8, &
172 : 0.96504E-08_r8, 0.18494E-07_r8, 0.34296E-07_r8, 0.61112E-07_r8, 0.10738E-06_r8, &
173 : 0.18747E-06_r8, 0.32054E-06_r8, 0.52872E-06_r8, 0.83634E-06_r8, 0.12723E-05_r8, &
174 : 0.18748E-05_r8, 0.26362E-05_r8, 0.36986E-05_r8, 0.52079E-05_r8, 0.72579E-05_r8, &
175 : 0.98614E-05_r8, 0.12775E-04_r8, 0.15295E-04_r8, 0.18072E-04_r8, 0.20959E-04_r8, &
176 : 0.19208E-04_r8, 0.16285E-04_r8, 0.13628E-04_r8, 0.11784E-04_r8, 0.11085E-04_r8, &
177 : 0.11916E-04_r8, 0.14771E-04_r8, 0.20471E-04_r8, 0.29426E-04_r8, 0.42992E-04_r8, &
178 : 0.62609E-04_r8, 0.90224E-04_r8, 0.12870E-03_r8, 0.18281E-03_r8, 0.26029E-03_r8, &
179 : 0.37224E-03_r8, 0.53254E-03_r8, 0.75697E-03_r8, 0.10623E-02_r8, 0.14660E-02_r8, &
180 : 0.19856E-02_r8, 0.26393E-02_r8, 0.34473E-02_r8, 0.44327E-02_r8, 0.56254E-02_r8, &
181 : 0.70672E-02_r8, 0.88174E-02_r8, 0.10960E-01_r8, 0.13613E-01_r8, 0.16934E-01_r8, &
182 : 0.21137E-01_r8, 0.26501E-01_r8, 0.33388E-01_r8, 0.42263E-01_r8, 0.53716E-01_r8, &
183 : 0.68491E-01_r8, 0.87521E-01_r8, 0.11196E+00_r8, 0.14320E+00_r8, 0.18295E+00_r8, &
184 : 0.23321E+00_r8, 0.29631E+00_r8/
185 :
186 : logical :: ionvels_read_from_file = .false.
187 :
188 : contains
189 :
190 : !==============================================================================
191 :
192 1536 : subroutine iondrag_register
193 : !-----------------------------------------------------------------------
194 : ! Register E and B fields.
195 : !
196 : ! Register iondrag variables with physics buffer:
197 : !
198 : ! Hall and Pedersen conductivities
199 : !
200 : ! pcols dimension and lchnk assumed here
201 : !
202 : !-----------------------------------------------------------------------
203 : use exbdrift, only: exbdrift_register
204 : use physics_buffer, only: pbuf_add_field, dtype_r8
205 :
206 : ! E and B fields
207 1536 : call exbdrift_register()
208 :
209 1536 : if ( waccmx_is("ionosphere") ) then
210 : ! Pedersen Conductivity and Hall Conductivity
211 0 : call pbuf_add_field('PedConduct', 'physpkg', dtype_r8, (/pcols,pver/), PedConduct_idx )
212 0 : call pbuf_add_field('HallConduct', 'physpkg', dtype_r8, (/pcols,pver/), HallConduct_idx)
213 : end if
214 :
215 1536 : end subroutine iondrag_register
216 :
217 : !================================================================================================
218 :
219 1536 : subroutine iondrag_readnl(nlfile)
220 :
221 1536 : use namelist_utils, only: find_group_name
222 : use units, only: getunit, freeunit
223 : use spmd_utils, only: mpicom, masterprocid, mpicom, mpi_character, mpi_logical, mpi_real8
224 :
225 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
226 :
227 : ! Local variables
228 : integer :: unitn, ierr
229 : character(len=*), parameter :: subname = 'iondrag_readnl'
230 :
231 : namelist /iondrag_nl/ efield_lflux_file, efield_hflux_file, efield_potential_max, empirical_ion_velocities
232 :
233 1536 : if (masterproc) then
234 2 : unitn = getunit()
235 2 : open( unitn, file=trim(nlfile), status='old' )
236 2 : call find_group_name(unitn, 'iondrag_nl', status=ierr)
237 2 : if (ierr == 0) then
238 2 : read(unitn, iondrag_nl, iostat=ierr)
239 2 : if (ierr /= 0) then
240 0 : call endrun(subname // ':: ERROR reading namelist')
241 : end if
242 : end if
243 2 : close(unitn)
244 2 : call freeunit(unitn)
245 :
246 : end if
247 :
248 1536 : call mpi_bcast (efield_lflux_file, len(efield_lflux_file), mpi_character, masterprocid, mpicom, ierr)
249 1536 : call mpi_bcast (efield_hflux_file, len(efield_hflux_file), mpi_character, masterprocid, mpicom, ierr)
250 1536 : call mpi_bcast (efield_potential_max, 1, mpi_real8, masterprocid, mpicom, ierr)
251 1536 : call mpi_bcast (empirical_ion_velocities, 1, mpi_logical, masterprocid, mpicom, ierr)
252 :
253 1536 : if (masterproc) then
254 2 : write(iulog,*) 'iondrag options:'
255 2 : write(iulog,*) ' efield_lflux_file = '//trim(efield_lflux_file)
256 2 : write(iulog,*) ' efield_hflux_file = '//trim(efield_hflux_file)
257 2 : write(iulog,*) ' efield_potential_max = ',efield_potential_max
258 2 : write(iulog,*) ' empirical_ion_velocities = ',empirical_ion_velocities
259 : end if
260 :
261 1536 : end subroutine iondrag_readnl
262 :
263 : !================================================================================================
264 :
265 1536 : subroutine iondrag_init( pref_mid )
266 : use constituents, only: cnst_get_ind
267 : use short_lived_species, only: slvd_index
268 :
269 : !-------------------------------------------------------------------------------
270 : ! Iondrag initialization, called from inti.F90.
271 : !-------------------------------------------------------------------------------
272 :
273 :
274 : !-------------------------------------------------------------------------------
275 : ! dummy arguments
276 : !-------------------------------------------------------------------------------
277 : real(r8), intent(in) :: pref_mid(pver)
278 :
279 : integer :: k, err
280 :
281 : integer :: cnst_ids(7)
282 :
283 1536 : doiodrg = .false.
284 1536 : do_waccm_ions = .false.
285 :
286 : !-------------------------------------------------------------------------------
287 : ! find lower bnd for iondrag
288 : !-------------------------------------------------------------------------------
289 1536 : if( pref_mid(1) < 0.1_r8 ) then
290 109056 : do k = 1, pver
291 109056 : if (pref_mid(k) < 50._r8) nbot_lev = k
292 : end do
293 : end if
294 :
295 1536 : if (nbot_lev > 0) then
296 1536 : doiodrg = .true.
297 : endif
298 :
299 1536 : if ( .not. doiodrg .and. masterproc ) then
300 0 : write(iulog,*) ' '
301 0 : write(iulog,*) 'iondrag_init: Does not have waccm level. Ion drag does not apply. '
302 0 : write(iulog,*) ' '
303 0 : return
304 : endif
305 :
306 1536 : if( masterproc ) then
307 2 : write(iulog,*) ' '
308 2 : write(iulog,*) 'iondrag_init: nbot_lev,press = ',nbot_lev,pref_mid(nbot_lev)
309 2 : write(iulog,*) ' '
310 : end if
311 :
312 1536 : call cnst_get_ind( 'e', id_elec, abort=.false. )
313 1536 : if (id_elec < 0) then
314 1536 : id_elec = slvd_index( 'e' )
315 : endif
316 1536 : call cnst_get_ind( 'Op', id_op, abort=.false. )
317 1536 : if (id_op < 0) then
318 1536 : id_op = slvd_index( 'Op' )
319 1536 : if (id_op > 0) then
320 1536 : op_slvd = .true.
321 : endif
322 : else
323 0 : op_slvd = .false.
324 : endif
325 1536 : call cnst_get_ind( 'O2p', id_o2p, abort=.false. )
326 1536 : if (id_o2p < 0) then
327 1536 : id_o2p = slvd_index( 'O2p' )
328 1536 : if (id_o2p > 0) then
329 1536 : o2p_slvd = .true.
330 : endif
331 : else
332 0 : o2p_slvd = .false.
333 : endif
334 1536 : call cnst_get_ind( 'NOp', id_nop, abort=.false. )
335 1536 : if (id_nop < 0) then
336 1536 : id_nop = slvd_index( 'NOp' )
337 1536 : if (id_nop > 0) then
338 1536 : nop_slvd = .true.
339 : endif
340 : else
341 0 : nop_slvd = .false.
342 : endif
343 1536 : call cnst_get_ind( 'O', id_xo1, abort=.false. )
344 1536 : if (id_xo1 < 0) then
345 0 : id_xo1 = slvd_index( 'O' )
346 0 : if (id_xo1 > 0) then
347 0 : xo1_slvd = .true.
348 : endif
349 : else
350 1536 : xo1_slvd = .false.
351 : endif
352 1536 : call cnst_get_ind( 'O2', id_xo2, abort=.false. )
353 1536 : if (id_xo2 < 0) then
354 0 : id_xo2 = slvd_index( 'O2' )
355 0 : if (id_xo2 > 0) then
356 0 : xo2_slvd = .true.
357 : endif
358 : else
359 1536 : xo2_slvd = .false.
360 : endif
361 1536 : call cnst_get_ind( 'N', id_n, abort=.false. )
362 1536 : if (id_n < 0) then
363 0 : id_n = slvd_index( 'N' )
364 : endif
365 :
366 12288 : cnst_ids = (/ id_elec, id_op, id_o2p, id_nop, id_xo1, id_xo2, id_n /)
367 :
368 12288 : if ( all( cnst_ids > 0 ) ) then
369 1536 : do_waccm_ions = .true.
370 : endif
371 :
372 1536 : if ( do_waccm_ions ) then
373 1536 : call ions_init
374 : else
375 0 : call ghg_init(pref_mid)
376 : endif
377 :
378 1536 : if (.not.doiodrg) return
379 :
380 : ! Add to masterfield list
381 3072 : call addfld('UIONTEND',(/ 'lev' /),'A','M/S2','u-tendency due to ion drag')
382 3072 : call addfld('VIONTEND',(/ 'lev' /),'A','M/S2','v-tendency due to ion drag')
383 :
384 : !
385 : ! Indices to 3d ion drifts :
386 1536 : ui_idx = pbuf_get_index('UI')
387 1536 : vi_idx = pbuf_get_index('VI')
388 1536 : wi_idx = pbuf_get_index('WI')
389 :
390 1536 : indxTe = pbuf_get_index( 'TElec',errcode=err )
391 1536 : indxTi = pbuf_get_index( 'TIon',errcode=err )
392 :
393 : end subroutine iondrag_init
394 :
395 : !================================================================================================
396 1536 : subroutine ions_init
397 :
398 : use efield, only: efield_init
399 : use exbdrift, only: exbdrift_init
400 : use mo_chem_utls, only: get_spc_ndx
401 : use chem_mods, only: adv_mass
402 : use mo_chem_utls, only: get_inv_ndx
403 : use chem_mods, only: fix_mass
404 : use phys_control, only: phys_getopts
405 :
406 : !-------------------------------------------------------------------------------
407 : ! local variables
408 : !-------------------------------------------------------------------------------
409 : integer :: id
410 :
411 : logical :: history_waccm
412 :
413 1536 : call phys_getopts(history_waccm_out=history_waccm)
414 :
415 : !-------------------------------------------------------------------------------
416 : ! initialize related packages: electric field
417 : !-------------------------------------------------------------------------------
418 :
419 1536 : call efield_init (efield_lflux_file, efield_hflux_file, efield_potential_max)
420 1536 : call exbdrift_init( empirical_ion_velocities )
421 :
422 1536 : id = get_spc_ndx('Op')
423 1536 : rmass_op = adv_mass(id)
424 1536 : id = get_spc_ndx('O2p')
425 1536 : rmass_o2p = adv_mass(id)
426 1536 : id = get_spc_ndx('NOp')
427 1536 : rmass_nop = adv_mass(id)
428 1536 : id = get_spc_ndx('O')
429 1536 : rmass_o1 = adv_mass(id)
430 1536 : id = get_spc_ndx('O2')
431 1536 : rmass_o2 = adv_mass(id)
432 1536 : id = get_inv_ndx('N2')
433 1536 : rmass_n2 = fix_mass(id)
434 :
435 1536 : rmi_o1 = 1._r8/rmass_o1
436 1536 : rmi_o2 = 1._r8/rmass_o2
437 1536 : rmi_n2 = 1._r8/rmass_n2
438 1536 : rmi_op = 1._r8/rmass_op
439 1536 : rmi_o2p = 1._r8/rmass_o2p
440 1536 : rmi_nop = 1._r8/rmass_nop
441 1536 : rmi_op_kg = 1._r8/(rmass_op *amu)
442 1536 : rmi_o2p_kg = 1._r8/(rmass_o2p*amu)
443 1536 : rmi_nop_kg = 1._r8/(rmass_nop*amu)
444 :
445 : !-------------------------------------------------------------------------------
446 : ! Set up fields to history files.
447 : !-------------------------------------------------------------------------------
448 1536 : call addfld('BNORTH', horiz_only, 'I', 'nT', 'northward component of magnetic field (nT)')
449 1536 : call addfld('BEAST' , horiz_only, 'I', 'nT', 'eastward component of magnetic field (nT)')
450 1536 : call addfld('BDOWN' , horiz_only, 'I', 'nT', 'downward component of magnetic field (nT)')
451 1536 : call addfld('BMAG' , horiz_only, 'I', 'nT', 'magnetic field magnitude (nT)' )
452 :
453 3072 : call addfld('QIONSUM', (/ 'lev' /), 'I','S-1' ,'Ion prod sum')
454 3072 : call addfld('ELECDEN', (/ 'lev' /), 'I','CM-3','NE (ion sum)')
455 3072 : call addfld('SIGMAPED', (/ 'lev' /), 'I', 'siemens/m', 'Pederson conductivity' )
456 3072 : call addfld('SIGMAHAL', (/ 'lev' /), 'I', 'siemens/m', 'Hall conductivity' )
457 3072 : call addfld('LAMDA1' , (/ 'lev' /), 'I' ,'S-1','LAMDA PED')
458 3072 : call addfld('LAMDA2' , (/ 'lev' /), 'I' ,'S-1','LAMDA HALL')
459 :
460 3072 : call addfld('LXX', (/ 'lev' /), 'I','S-1','LXX')
461 3072 : call addfld('LYY', (/ 'lev' /), 'I','S-1','LYY')
462 3072 : call addfld('LXY', (/ 'lev' /), 'I','S-1','LXY')
463 3072 : call addfld('LYX', (/ 'lev' /), 'I','S-1','LYX')
464 :
465 : !
466 : ! Joule heating, and tn before and after joule heating tendencies are applied:
467 : !
468 3072 : call addfld( 'QJOULE', (/ 'lev' /), 'I', 'K/s' , 'Joule Heat' ) ! joule heating
469 1536 : if (history_waccm) then
470 1536 : call add_default( 'QJOULE ', 1, ' ' ) ! joule heating (K/s)
471 : end if
472 : !
473 : ! 3d drifts from either edynamo or exbdrift.
474 : !
475 1536 : if (empirical_ion_velocities) then
476 3072 : call addfld('UI',(/ 'lev' /),'I','m/s', 'UI Zonal empirical ExB drift from exbdrift')
477 3072 : call addfld('VI',(/ 'lev' /),'I','m/s', 'VI Meridional empirical ExB drift from exbdrift')
478 3072 : call addfld('WI',(/ 'lev' /),'I','m/s', 'WI Vertical empirical ExB drift from exbdrift')
479 : endif
480 :
481 1536 : end subroutine ions_init
482 :
483 : !========================================================================
484 :
485 0 : subroutine ghg_init (pref_mid)
486 :
487 : !
488 : ! initialization for ion drag calculation
489 : !
490 :
491 : !------------------Input arguments---------------------------------------
492 :
493 : real(r8), intent(in) :: pref_mid(pver) ! model ref pressure at midpoint
494 :
495 : !-----------------local workspace---------------------------------------
496 : integer k
497 : integer kinv
498 :
499 : real(r8) rpsh ! ref pressure scale height
500 :
501 : real(r8), parameter :: preftgcm = 5.e-5_r8 ! TIME GCM reference pressure (Pa)
502 :
503 : !------------------------------------------------------------------------
504 :
505 : ! With the defualt values of nbot_lev and ntop_lev, ion drag calcualtion are NOT carried out
506 0 : nbot_lev=0
507 0 : ntop_lev=1
508 :
509 0 : do k = 1, pver
510 0 : rpsh=log(1e5_r8/pref_mid(k))
511 0 : if (rpsh .gt. 14._r8) nbot_lev = k
512 : end do
513 0 : if (nbot_lev .gt. ntop_lev) doiodrg=.true.
514 0 : if (masterproc) then
515 0 : write(iulog,fmt='(a15)') 'From IONDRAGI:'
516 0 : write(iulog,fmt='(1a12,1i10)') 'NTOP_LEV =',ntop_lev
517 0 : write(iulog,fmt='(1a12,1i10)') 'NBOT_LEV =',nbot_lev
518 0 : write(iulog,*) 'IONDRAG flag is',doiodrg
519 : endif
520 0 : if (.not.doiodrg) return
521 :
522 : ! obtain TIME/GCM pressure scale height
523 0 : pshtiod(1)=-17._r8
524 0 : do k=2,plevtiod
525 0 : pshtiod(k)=pshtiod(k-1)+0.25_r8
526 : enddo
527 :
528 : ! map TIME-psh into CCM-psh
529 0 : pshtiod=pshtiod-log(preftgcm/1E5_r8)
530 :
531 : ! CCM psh
532 : ! note that vertical indexing is inverted with respect to CCM standard
533 0 : do k=1,pver
534 0 : kinv=pver-k+1
535 0 : pshccm(kinv)=log(1e5_r8/pref_mid(k))
536 : enddo
537 :
538 : ! vertical interpolation
539 0 : write(iulog,*) ' '
540 0 : write(iulog,*) 'iondragi: before lininterp for alamxx'
541 0 : write(iulog,*) ' nlatin,nlatout =',plevtiod,pver
542 0 : write(iulog,*) ' yin'
543 0 : write(iulog,'(1p,5g15.8)') pshtiod
544 0 : write(iulog,*) ' yout'
545 0 : write(iulog,'(1p,5g15.8)') pshccm
546 0 : write(iulog,*) ' '
547 :
548 0 : call lininterp (alamxx ,pshtiod,plevtiod, alamxxi ,pshccm,pver)
549 :
550 0 : call lininterp (alamxy ,pshtiod,plevtiod, alamxyi ,pshccm,pver)
551 :
552 : ! invert indeces back to CCM convention
553 0 : alamxxi(1:pver)=alamxxi(pver:1:-1)
554 0 : alamxyi(1:pver)=alamxyi(pver:1:-1)
555 :
556 0 : return
557 :
558 1536 : end subroutine ghg_init
559 :
560 : !================================================================================================
561 :
562 16128 : subroutine iondrag_timestep_init
563 : use efield, only: get_efield
564 :
565 16128 : if (do_waccm_ions) then
566 16128 : if (empirical_ion_velocities .or. (is_first_step().and..not.ionvels_read_from_file)) then
567 : ! Compute the electric field
568 16128 : call t_startf ('efield')
569 16128 : call get_efield
570 16128 : call t_stopf ('efield')
571 : endif
572 : endif
573 :
574 16128 : end subroutine iondrag_timestep_init
575 :
576 : !================================================================================================
577 15978240 : subroutine iondrag_calc_ions( lchnk, ncol, state, ptend, pbuf, delt )
578 : !-------------------------------------------------------------------------------
579 : ! Calculate ion drag tensors lxx,lyy,lxy,lyx.
580 : ! Also calculate Pedersen and Hall conductivities.
581 : ! This is called from tphysac.
582 : !-------------------------------------------------------------------------------
583 :
584 : use mo_apex,only: & ! (pcols,begchunk:endchunk)
585 : bnorth, & ! northward component of magnetic field (nT)
586 : beast, & ! eastward component of magnetic field (nT)
587 16128 : bdown, & ! downward component of magnetic field (nT)
588 : bmag ! magnetic field magnitude (nT)
589 : use physconst, only: avogad, boltz
590 : use chemistry, only: imozart
591 : use mo_mean_mass, only: set_mean_mass
592 : use exbdrift, only: exbdrift_ion_vels
593 : use short_lived_species, only: slvd_pbf_ndx => pbf_idx
594 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
595 :
596 : !-------------------------------------------------------------------------------
597 : ! dummy arguments
598 : !-------------------------------------------------------------------------------
599 : integer,intent(in) :: lchnk ! current chunk index
600 : integer,intent(in) :: ncol ! number of atmospheric columns
601 : real(r8), intent(in) :: delt ! time step (s)
602 : type(physics_state), intent(in), target :: state ! Physics state variables
603 : type(physics_ptend), intent(out) :: ptend ! Physics tendencies
604 :
605 : type(physics_buffer_desc), pointer :: pbuf(:)
606 :
607 : !-------------------------------------------------------------------------------
608 : ! Local variables
609 : !-------------------------------------------------------------------------------
610 : integer :: i,k ! loop indices
611 :
612 : real(r8) :: sqrt_te ! sqrt(te)
613 : real(r8) :: sqrt_tnti ! sqrt(tnti)
614 : real(r8) :: wrk
615 : real(r8),parameter :: dipmin = 0.17_r8 ! minimum dip angle (tuneable)
616 : real(r8),parameter :: emass = 9.1093819e-31_r8 ! electron mass (kg)
617 : real(r8),parameter :: qe = 1.6021765e-19_r8 ! electronic charge (coulombs)
618 : real(r8),parameter :: colfac = 1.5_r8 ! collision factor (tuneable)
619 : real(r8),parameter :: boltzmann = 1.e7_r8 * boltz
620 : real(r8),parameter :: avo = avogad*1.e-3_r8 ! (molecules/mole)
621 : ! real(r8),parameter :: rmass_op = 15.9989_r8 ! mass of O+
622 : ! real(r8),parameter :: rmass_o2p = 31.9983_r8 ! mass of O2+
623 : ! real(r8),parameter :: rmass_nop = 30.0056_r8 ! mass of NO+
624 : ! real(r8),parameter :: rmass_o1 = 16._r8 ! mass of O
625 : ! real(r8),parameter :: rmass_o2 = 32._r8 ! mass of O2
626 : ! real(r8),parameter :: rmass_n2 = 28._r8 ! mass of N2
627 :
628 : !-------------------------------------------------------------------------------
629 : ! Inverted masses (for multiply in loops rather than divide):
630 : !-------------------------------------------------------------------------------
631 : ! real(r8),parameter :: rmi_o1 = 1._r8/rmass_o1
632 : ! real(r8),parameter :: rmi_o2 = 1._r8/rmass_o2
633 : ! real(r8),parameter :: rmi_n2 = 1._r8/rmass_n2
634 : ! real(r8),parameter :: rmi_op = 1._r8/rmass_op
635 : ! real(r8),parameter :: rmi_o2p = 1._r8/rmass_o2p
636 : ! real(r8),parameter :: rmi_nop = 1._r8/rmass_nop
637 : ! real(r8),parameter :: rmi_op_kg = 1._r8/(rmass_op *amu)
638 : ! real(r8),parameter :: rmi_o2p_kg = 1._r8/(rmass_o2p*amu)
639 : ! real(r8),parameter :: rmi_nop_kg = 1._r8/(rmass_nop*amu)
640 :
641 : real(r8), target :: tn(pcols,pver) ! neutral gas temperature (deg K)
642 : real(r8) :: xo2 (pcols,pver) ! O2 (mmr)
643 : real(r8) :: xo1 (pcols,pver) ! O (mmr)
644 : real(r8) :: xn2 (pcols,pver) ! N2 (mmr)
645 : real(r8) :: o2p (pcols,pver) ! O2+ (mmr)
646 : real(r8) :: op (pcols,pver) ! O+ (mmr)
647 : real(r8) :: nop (pcols,pver) ! NO+ (mmr)
648 : real(r8) :: barm (pcols,pver) ! mean molecular weight (g/mole)
649 : real(r8) :: xnmbar (pcols,pver) ! for unit conversion to volume density
650 : real(r8) :: tnti (pcols,pver) ! average of tn and ti
651 : real(r8) :: o2_cm3 (pcols,pver) ! o2 volume density (cm-3)
652 : real(r8) :: o1_cm3 (pcols,pver) ! o volume density (cm-3)
653 : real(r8) :: n2_cm3 (pcols,pver) ! n2 volume density (cm-3)
654 : real(r8) :: o2p_cm3 (pcols,pver) ! O2+ (cm-3)
655 : real(r8) :: op_cm3 (pcols,pver) ! O+ (cm-3)
656 : real(r8) :: nop_cm3 (pcols,pver) ! NO+ (cm-3)
657 : real(r8) :: ne_sigmas(pcols,pver) ! electron density for conductivities
658 : real(r8) :: lamda1 (pcols,pver) ! sigped*b**2/rho
659 : real(r8) :: lamda2 (pcols,pver) ! sighal*b**2/rho
660 : real(r8) :: lxxnorot(pcols,pver) ! XX before rotation
661 : real(r8) :: lyynorot(pcols,pver) ! YY before rotation
662 : real(r8) :: lxynorot(pcols,pver) ! XY before rotation
663 :
664 : !-------------------------------------------------------------------------------
665 : ! Ion-neutral momentum transfer collision frequencies.
666 : ! rnu_xxx_xx = ratio of collision to gyro-frequences for O2+, O+, NO+.
667 : !-------------------------------------------------------------------------------
668 : real(r8) :: rnu_o2p_o2(pcols,pver) ! O2+ ~ O2 collision freq (resonant, T dependent)
669 : real(r8) :: rnu_op_o2 (pcols,pver) ! O+ ~ O2 collision freq (non-resonant)
670 : real(r8) :: rnu_nop_o2(pcols,pver) ! NO+ ~ O2 collision freq (non-resonant)
671 : real(r8) :: rnu_o2p_o (pcols,pver) ! O2+ ~ O collision freq (non-resonant)
672 : real(r8) :: rnu_op_o (pcols,pver) ! O+ ~ O collision freq (resonant, T dependent)
673 : real(r8) :: rnu_nop_o (pcols,pver) ! NO+ ~ O collision freq (non-resonant)
674 : real(r8) :: rnu_o2p_n2(pcols,pver) ! O2+ ~ N2 collision freq (non-resonant)
675 : real(r8) :: rnu_op_n2 (pcols,pver) ! O+ ~ N2 collision freq (non-resonant)
676 : real(r8) :: rnu_nop_n2(pcols,pver) ! NO+ ~ N2 collision freq (non-resonant)
677 : real(r8) :: rnu_o2p (pcols,pver) ! [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
678 : real(r8) :: rnu_op (pcols,pver) ! [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
679 : real(r8) :: rnu_nop (pcols,pver) ! [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
680 : real(r8) :: rnu_ne (pcols,pver) ! electron ~ neutral collision frequency (s-1)
681 :
682 : real(r8) :: press (pcols) ! pressure at interface levels (dyne/cm^2)
683 : real(r8) :: qe_fac (pcols) ! unit conversion factor for conductivities
684 : real(r8) :: dipmag (pcols) ! magnetic dip angle
685 : real(r8) :: decmag (pcols) ! magnetic declination
686 : real(r8) :: btesla (pcols) ! magnetic field (teslas)
687 : real(r8) :: sindip (pcols) ! sin(dipmag)
688 : real(r8) :: sin2dip (pcols) ! sindip^2
689 : real(r8) :: sindec (pcols) ! sin(decmag)
690 : real(r8) :: cosdec (pcols) ! cos(decmag)
691 : real(r8) :: sin2dec (pcols) ! sindec^2
692 : real(r8) :: cos2dec (pcols) ! cosdec^2
693 : real(r8) :: omega_o2p (pcols) ! angular gyrofrequency for o2+ (s-1)
694 : real(r8) :: omega_op (pcols) ! angular gyrofrequency for o+ (s-1)
695 : real(r8) :: omega_nop (pcols) ! angular gyrofrequency for no+ (s-1)
696 : real(r8) :: omega_e (pcols) ! electron angular gyrofrequency (s-1)
697 : real(r8) :: omega_o2p_inv(pcols) ! inverse of o2+ gyrofrequency
698 : real(r8) :: omega_op_inv (pcols) ! inverse of o+ gyrofrequency
699 : real(r8) :: omega_nop_inv(pcols) ! inverse of no+ gyrofrequency
700 : real(r8) :: omega_e_inv (pcols) ! inverse of electron gyrofrequency
701 :
702 : !-------------------------------------------------------------------------------
703 : ! Ion drag coefficients output:
704 : !-------------------------------------------------------------------------------
705 : real(r8) :: lxx(pcols,pver) ! lambda XX coefficients (s-1)
706 : real(r8) :: lyy(pcols,pver) ! lambda YY coefficients (s-1)
707 : real(r8) :: lxy(pcols,pver) ! lambda XY coefficients (s-1)
708 : real(r8) :: lyx(pcols,pver) ! lambda YX coefficients (s-1)
709 :
710 : !-------------------------------------------------------------------------------
711 : ! Conductivities output:
712 : !-------------------------------------------------------------------------------
713 : real(r8) :: sigma_ped (pcols,pver) ! pedersen conductivity (siemens/m)
714 : real(r8) :: sigma_hall(pcols,pver) ! hall conductivity (siemens/m)
715 :
716 : real(r8) :: qout(pcols,pver) ! temp for outfld
717 :
718 72960 : real(r8), dimension(:,:), pointer :: q_xo1, q_xo2, q_o2p, q_op, q_nop
719 :
720 72960 : real(r8), dimension(:,:), pointer :: tE ! electron temperature in pbuf (K)
721 72960 : real(r8), dimension(:,:), pointer :: tI ! ion temperature in pbuf (K)
722 :
723 72960 : if (.not.doiodrg) return
724 :
725 72960 : call outfld ('BNORTH', bnorth(:,lchnk), pcols, lchnk )
726 72960 : call outfld ('BEAST' , beast(:,lchnk), pcols, lchnk )
727 72960 : call outfld ('BDOWN' , bdown(:,lchnk), pcols, lchnk )
728 72960 : call outfld ('BMAG' , bmag(:,lchnk), pcols, lchnk )
729 :
730 72960 : if ( xo1_slvd ) then
731 0 : call pbuf_get_field(pbuf, slvd_pbf_ndx, q_xo1, start=(/1,1,id_xo1/), kount=(/pcols,pver,1/) )
732 : else
733 72960 : q_xo1 => state%q(:,:,id_xo1)
734 : endif
735 72960 : if ( xo2_slvd ) then
736 0 : call pbuf_get_field(pbuf, slvd_pbf_ndx, q_xo2, start=(/1,1,id_xo2/), kount=(/pcols,pver,1/) )
737 : else
738 72960 : q_xo2 => state%q(:,:,id_xo2)
739 : endif
740 72960 : if ( o2p_slvd ) then
741 291840 : call pbuf_get_field(pbuf, slvd_pbf_ndx, q_o2p, start=(/1,1,id_o2p/), kount=(/pcols,pver,1/) )
742 : else
743 0 : q_o2p => state%q(:,:,id_o2p)
744 : endif
745 72960 : if ( op_slvd ) then
746 291840 : call pbuf_get_field(pbuf, slvd_pbf_ndx, q_op, start=(/1,1,id_op/), kount=(/pcols,pver,1/) )
747 : else
748 0 : q_op => state%q(:,:,id_op)
749 : endif
750 72960 : if ( nop_slvd ) then
751 291840 : call pbuf_get_field(pbuf, slvd_pbf_ndx, q_nop, start=(/1,1,id_nop/), kount=(/pcols,pver,1/) )
752 : else
753 0 : q_nop => state%q(:,:,id_nop)
754 : endif
755 :
756 : !-------------------------------------------------------------------------------
757 : ! Define local tn and major species from state (mmr):
758 : !-------------------------------------------------------------------------------
759 5180160 : do k = 1,pver
760 78723840 : do i = 1,ncol
761 73543680 : tn (i,k) = state%t(i,k)
762 73543680 : xo2(i,k) = q_xo2(i,k) ! o2 (mmr)
763 73543680 : xo1(i,k) = q_xo1(i,k) ! o (mmr)
764 73543680 : xn2(i,k) = 1._r8 - (xo2(i,k) + xo1(i,k)) ! n2 (mmr)
765 73543680 : xn2(i,k) = max( 1.e-20_r8,xn2(i,k) )
766 73543680 : o2p(i,k) = q_o2p(i,k) ! o2+ (mmr)
767 73543680 : op (i,k) = q_op(i,k) ! o+ (mmr)
768 78650880 : nop(i,k) = q_nop(i,k) ! no+ (mmr)
769 : end do
770 : end do
771 :
772 : !--------------------------------------------------------------------------------------------------
773 : ! For WACCM-X, grab electron (TE) and ion (TI) temperatures from pbuf from ionosphere module
774 : !--------------------------------------------------------------------------------------------------
775 72960 : if ( indxTe>0 .and. indxTi>0 ) then
776 0 : call pbuf_get_field(pbuf, indxTe, tE)
777 0 : call pbuf_get_field(pbuf, indxTi, tI)
778 : else
779 72960 : tE => tn
780 72960 : tI => tn
781 : endif
782 :
783 78723840 : qout(:ncol,:) = o2p(:ncol,:) + op(:ncol,:) + nop(:ncol,:)
784 72960 : call outfld ('QIONSUM ', qout, pcols, lchnk)
785 :
786 : !-------------------------------------------------------------------------------
787 : ! calculate empirical ExB drift velocities if not using edynamo drifts
788 : ! (if edynamo calculated the drifts, then they are already in pbuf%ui, etc.)
789 : !-------------------------------------------------------------------------------
790 72960 : if (empirical_ion_velocities .or. (is_first_step().and..not.ionvels_read_from_file)) then
791 72960 : call t_startf ( 'exbdrift_ion_vels' )
792 72960 : call exbdrift_ion_vels( lchnk, ncol, pbuf)
793 72960 : call t_stopf ( 'exbdrift_ion_vels' )
794 : endif
795 :
796 : !-------------------------------------------------------------------------------
797 :
798 1123584 : do i = 1,ncol
799 1050624 : btesla(i) = bmag(i,lchnk)*1.e-9_r8 ! nT to teslas (see bmag in apex module)
800 : !-------------------------------------------------------------------------------
801 : ! Angular gyrofrequency of O+, O2+ and NO+ (s-1):
802 : !-------------------------------------------------------------------------------
803 1050624 : omega_op (i) = qe*btesla(i)*rmi_op_kg
804 1050624 : omega_o2p(i) = qe*btesla(i)*rmi_o2p_kg
805 1050624 : omega_nop(i) = qe*btesla(i)*rmi_nop_kg
806 : !-------------------------------------------------------------------------------
807 : ! Electron angular gyrofrequency (s-1):
808 : !-------------------------------------------------------------------------------
809 1050624 : omega_e(i) = qe*btesla(i)/emass
810 : !-------------------------------------------------------------------------------
811 : ! Invert now, so we can multiply rather than divide in loops below:
812 : !-------------------------------------------------------------------------------
813 1050624 : omega_op_inv (i) = 1._r8/omega_op(i)
814 1050624 : omega_o2p_inv(i) = 1._r8/omega_o2p(i)
815 1050624 : omega_nop_inv(i) = 1._r8/omega_nop(i)
816 1050624 : omega_e_inv(i) = 1._r8/omega_e(i)
817 : !-------------------------------------------------------------------------------
818 : ! Magnetic field geometry (used below in rotation of lambdas):
819 : !-------------------------------------------------------------------------------
820 1050624 : dipmag(i) = atan( bdown(i,lchnk)/sqrt(bnorth(i,lchnk)**2+beast(i,lchnk)**2) )
821 1050624 : decmag(i) = -atan2( beast(i,lchnk),bnorth(i,lchnk) )
822 1050624 : cosdec(i) = cos( decmag(i) )
823 1050624 : sindec(i) = sin( decmag(i) )
824 1050624 : if( abs(dipmag(i)) >= dipmin ) then
825 997139 : sindip(i) = sin(dipmag(i))
826 : else
827 53485 : if( dipmag(i) >= 0._r8 ) then
828 26505 : sindip(i) = sin( dipmin )
829 : else
830 26980 : sindip(i) = sin( -dipmin )
831 : end if
832 : end if
833 1050624 : sin2dip(i) = sindip(i)**2
834 1050624 : sin2dec(i) = sindec(i)**2
835 1123584 : cos2dec(i) = cosdec(i)**2
836 : end do
837 :
838 : ! write(iulog,"('iondrag: btesla=', /,(6e12.4))") btesla
839 : ! write(iulog,"('iondrag: bdown=', /,(6e12.4))") bdown(:,lchnk)
840 : ! write(iulog,"('iondrag: beast=', /,(6e12.4))") beast(:,lchnk)
841 : ! write(iulog,"('iondrag: bnorth=', /,(6e12.4))") bnorth(:,lchnk)
842 : ! write(iulog,"('iondrag: omega_o2p=',/,(6e12.4))") omega_o2p
843 : ! write(iulog,"('iondrag: omega_op=' ,/,(6e12.4))") omega_op
844 : ! write(iulog,"('iondrag: omega_nop=',/,(6e12.4))") omega_nop
845 :
846 : !-------------------------------------------------------------------------------
847 : ! Ion-neutral momentum transfer collision frequency coefficients:
848 : !-------------------------------------------------------------------------------
849 5180160 : do k = 1,pver
850 78723840 : do i = 1,ncol
851 73543680 : tnti(i,k) = 0.5_r8*(tI(i,k) + tn(i,k)) ! ave of tn & ti
852 73543680 : sqrt_tnti = sqrt( tnti(i,k) )
853 73543680 : wrk = log10( tnti(i,k) )
854 : !-------------------------------------------------------------------------------
855 : ! Collision frequency coefficients with O2 (cm3/s):
856 : !-------------------------------------------------------------------------------
857 : rnu_o2p_o2(i,k) = 2.59e-11_r8*sqrt_tnti & ! O2+ ~ O2 (resonant)
858 73543680 : *(1._r8 - .073_r8*wrk)**2
859 73543680 : rnu_op_o2 (i,k) = 6.64e-10_r8 ! O+ ~ O2
860 73543680 : rnu_nop_o2(i,k) = 4.27e-10_r8 ! NO+ ~ O2
861 : !-------------------------------------------------------------------------------
862 : ! Collision frequency coefficients with O (cm3/s):
863 : !-------------------------------------------------------------------------------
864 73543680 : rnu_o2p_o(i,k) = 2.31e-10_r8 ! O2+ ~ O
865 : rnu_op_o (i,k) = 3.67e-11_r8*sqrt_tnti & ! O+ ~ O (resonant)
866 73543680 : *(1._r8 - .064_r8*wrk)**2*colfac
867 73543680 : rnu_nop_o(i,k) = 2.44e-10_r8 ! NO+ ~ O
868 : !-------------------------------------------------------------------------------
869 : ! Collision frequency coefficients with N2 (cm3/s):
870 : !-------------------------------------------------------------------------------
871 73543680 : rnu_o2p_n2(i,k) = 4.13e-10_r8 ! O2+ ~ N2
872 73543680 : rnu_op_n2 (i,k) = 6.82e-10_r8 ! O+ ~ N2
873 78650880 : rnu_nop_n2(i,k) = 4.34e-10_r8 ! NO+ ~ N2
874 : end do
875 : end do
876 :
877 : !-------------------------------------------------------------------------------
878 : ! Sub set_mean_mass (mo_mean_mass.F90) returns barm(ncol,pver) in g/mole,
879 : ! however, set_mean_mass sometimes returns zero in top(?) four values
880 : ! of the column, so barm is calculated here, see below.
881 : !
882 : ! call set_mean_mass(ncol, state%q(1,1,imozart), barm)
883 : !
884 : ! Major species and ion number densities (mmr to cm-3):
885 : !-------------------------------------------------------------------------------
886 :
887 72960 : call set_mean_mass( ncol, state%q(:,:,imozart:), barm )
888 :
889 5180160 : do k = 1,pver
890 78723840 : do i = 1,ncol
891 73543680 : press(i) = 10._r8*state%pmid(i,k) ! from Pa to dyne/cm^2
892 : ! barm(i,k) = 1._r8 / (xo2(i,k)*rmi_o2 + xo1(i,k)*rmi_o1 + xn2(i,k)*rmi_n2)
893 73543680 : xnmbar(i,k) = press(i)*barm(i,k)/(boltzmann*tn(i,k))
894 73543680 : o2_cm3(i,k) = xo2(i,k)*xnmbar(i,k)*rmi_o2 ! o2 (cm-3)
895 73543680 : o1_cm3(i,k) = xo1(i,k)*xnmbar(i,k)*rmi_o1 ! o (cm-3)
896 73543680 : n2_cm3(i,k) = xn2(i,k)*xnmbar(i,k)*rmi_n2 ! n2 (cm-3)
897 73543680 : o2p_cm3(i,k) = o2p(i,k)*xnmbar(i,k)*rmi_o2p ! o2+ (cm-3)
898 73543680 : op_cm3 (i,k) = op (i,k)*xnmbar(i,k)*rmi_op ! o+ (cm-3)
899 73543680 : nop_cm3(i,k) = nop(i,k)*xnmbar(i,k)*rmi_nop ! no+ (cm-3)
900 : !
901 : !----------------------------------------------------------------------------------
902 : ! Use sum of the 3 major ion number densities (as in tiegcm)
903 : !----------------------------------------------------------------------------------
904 : !
905 78650880 : ne_sigmas(i,k) = op_cm3(i,k) + o2p_cm3(i,k) + nop_cm3(i,k)
906 : end do
907 : end do
908 :
909 : !-------------------------------------------------------------------------------
910 : ! Multiply collision freq by neutral number density and sum for each ion:
911 : !
912 : ! rnu_o2p = [[o2p~o2]n(o2)+[o2p~o]n(o)+[o2p~n2]n(n2)]/w(o2p)
913 : ! rnu_op = [[op ~o2]n(o2)+[op ~o]n(o)+[op ~n2]n(n2)]/w(op )
914 : ! rnu_nop = [[nop~o2]n(o2)+[nop~o]n(o)+[nop~n2]n(n2)]/w(nop)
915 : !-------------------------------------------------------------------------------
916 5180160 : do k = 1,pver
917 78723840 : do i = 1,ncol
918 147087360 : rnu_o2p(i,k) = rnu_o2p_o2(i,k)*o2_cm3(i,k) &
919 : + rnu_o2p_o (i,k)*o1_cm3(i,k) &
920 147087360 : + rnu_o2p_n2(i,k)*n2_cm3(i,k)
921 : rnu_op (i,k) = rnu_op_o2 (i,k)*o2_cm3(i,k) &
922 : + rnu_op_o (i,k)*o1_cm3(i,k) &
923 73543680 : + rnu_op_n2 (i,k)*n2_cm3(i,k)
924 : rnu_nop(i,k) = rnu_nop_o2(i,k)*o2_cm3(i,k) &
925 : + rnu_nop_o (i,k)*o1_cm3(i,k) &
926 73543680 : + rnu_nop_n2(i,k)*n2_cm3(i,k)
927 : !-------------------------------------------------------------------------------
928 : ! Electron collision frequency (s-1):
929 : !-------------------------------------------------------------------------------
930 73543680 : sqrt_te = sqrt(tE(i,k))
931 : rnu_ne(i,k) = &
932 : 2.33e-11_r8*n2_cm3(i,k)*tE(i,k)*(1._r8 - 1.21e-4_r8*tE(i,k)) &
933 : + 1.82e-10_r8*o2_cm3(i,k)*sqrt_te*(1._r8 + 3.60e-2_r8*sqrt_te) &
934 78650880 : + 8.90e-11_r8*o1_cm3(i,k)*sqrt_te*(1._r8 + 5.70e-4_r8*tE(i,k))
935 : end do
936 : end do
937 :
938 : !-------------------------------------------------------------------------------
939 : ! Ratio of collision to gyro frequencies for o2+, o+, no+, ne:
940 : !-------------------------------------------------------------------------------
941 5180160 : do k = 1,pver
942 78723840 : do i = 1,ncol
943 73543680 : rnu_o2p(i,k) = rnu_o2p(i,k)*omega_o2p_inv(i)
944 73543680 : rnu_op (i,k) = rnu_op (i,k)*omega_op_inv (i)
945 73543680 : rnu_nop(i,k) = rnu_nop(i,k)*omega_nop_inv(i)
946 78650880 : rnu_ne (i,k) = rnu_ne (i,k)*omega_e_inv (i)
947 : end do
948 : end do
949 :
950 : !-------------------------------------------------------------------------------
951 : ! Calculate pedersen and Hall conductivities (siemens/m):
952 : !
953 : ! Qe_fac: 1.e6 to convert number densities from cm-3 to m-3:
954 : !-------------------------------------------------------------------------------
955 1123584 : qe_fac(:ncol) = qe*1.e6_r8/btesla(:ncol)
956 :
957 5180160 : do k = 1,pver
958 78723840 : do i = 1,ncol
959 :
960 : !-------------------------------------------------------------------------------
961 : ! Pedersen conductivity (siemens/m):
962 : !-------------------------------------------------------------------------------
963 147087360 : sigma_ped(i,k) = qe_fac(i) &
964 : *((op_cm3 (i,k)*rnu_op (i,k)/(1._r8 + rnu_op (i,k)**2)) &
965 : + (o2p_cm3 (i,k)*rnu_o2p(i,k)/(1._r8 + rnu_o2p(i,k)**2)) &
966 : + (nop_cm3 (i,k)*rnu_nop(i,k)/(1._r8 + rnu_nop(i,k)**2)) &
967 147087360 : + (ne_sigmas(i,k)*rnu_ne (i,k)/(1._r8 + rnu_ne (i,k)**2)))
968 : !-------------------------------------------------------------------------------
969 : ! Hall conductivity (siemens/m):
970 : !-------------------------------------------------------------------------------
971 : sigma_hall(i,k) = qe_fac(i) &
972 : *(ne_sigmas(i,k)/(1._r8 + rnu_ne (i,k)**2) &
973 : - op_cm3 (i,k)/(1._r8 + rnu_op (i,k)**2) &
974 : - o2p_cm3 (i,k)/(1._r8 + rnu_o2p(i,k)**2) &
975 78650880 : - nop_cm3 (i,k)/(1._r8 + rnu_nop(i,k)**2))
976 : end do
977 : end do
978 :
979 72960 : call outfld ('ELECDEN ',ne_sigmas ,pcols,lchnk)
980 72960 : call outfld ('SIGMAPED',sigma_ped ,pcols,lchnk)
981 72960 : call outfld ('SIGMAHAL',sigma_hall,pcols,lchnk)
982 :
983 : !--------------------------------------------------------------------------------------------
984 : ! Save conductivities in physics buffer using pointer for access in ionosphere module
985 : !--------------------------------------------------------------------------------------------
986 72960 : if ( waccmx_is('ionosphere') ) then
987 0 : call pbuf_set_field(pbuf, PedConduct_idx, sigma_ped(1:ncol,1:pver), start=(/1,1/), kount=(/ncol,pver/) )
988 0 : call pbuf_set_field(pbuf, HallConduct_idx, sigma_hall(1:ncol,1:pver), start=(/1,1/), kount=(/ncol,pver/) )
989 : endif
990 :
991 5180160 : do k = 1,pver
992 78723840 : do i = 1,ncol
993 73543680 : wrk = btesla(i)**2*avo*1.e-3_r8/xnmbar(i,k)
994 73543680 : lamda1(i,k) = sigma_ped(i,k)*wrk
995 78650880 : lamda2(i,k) = sigma_hall(i,k)*wrk
996 : end do
997 : end do
998 :
999 72960 : call outfld ('LAMDA1',lamda1,pcols,lchnk)
1000 72960 : call outfld ('LAMDA2',lamda2,pcols,lchnk)
1001 :
1002 5180160 : do k = 1,pver
1003 78723840 : do i = 1,ncol
1004 73543680 : lxxnorot(i,k) = lamda1(i,k)
1005 73543680 : lyynorot(i,k) = lamda1(i,k)*sin2dip(i)
1006 78650880 : lxynorot(i,k) = lamda2(i,k)*sindip(i)
1007 : end do
1008 : end do
1009 :
1010 : !-------------------------------------------------------------------------------
1011 : ! Rotate lambdas from local magnetic to geographic coordinates:
1012 : !-------------------------------------------------------------------------------
1013 5180160 : do k = 1,pver
1014 78723840 : do i = 1,ncol
1015 73543680 : lxx(i,k) = lxxnorot(i,k)*cos2dec(i) + lyynorot(i,k)*sin2dec(i)
1016 73543680 : lyy(i,k) = lyynorot(i,k)*cos2dec(i) + lxxnorot(i,k)*sin2dec(i)
1017 73543680 : wrk = (lyynorot(i,k) - lxxnorot(i,k))*sindec(i)*cosdec(i)
1018 73543680 : lyx(i,k) = lxynorot(i,k) - wrk
1019 78650880 : lxy(i,k) = lxynorot(i,k) + wrk
1020 : end do
1021 : end do
1022 :
1023 72960 : call outfld ('LXX ',lxx,pcols,lchnk)
1024 72960 : call outfld ('LYY ',lyy,pcols,lchnk)
1025 72960 : call outfld ('LXY ',lxy,pcols,lchnk)
1026 72960 : call outfld ('LYX ',lyx,pcols,lchnk)
1027 :
1028 72960 : call physics_ptend_init(ptend, state%psetcols, "ion drag", lu=.true., lv=.true., ls=.true.)
1029 :
1030 : !-------------------------------------------------------------------------------
1031 : ! Calculate ion drag tendencies and apply to neutral velocities:
1032 : !-------------------------------------------------------------------------------
1033 : call iondrag_tend( lchnk, ncol, state, ptend, pbuf, &
1034 72960 : lxx, lyy, lxy, lyx, delt )
1035 :
1036 : !-------------------------------------------------------------------------------
1037 : ! Calculate joule heating tendency and apply to temperature:
1038 : !-------------------------------------------------------------------------------
1039 : call jouleheat_tend( lchnk, ncol, state, ptend, pbuf, &
1040 72960 : lxx, lyy, lxy, lyx )
1041 :
1042 145920 : end subroutine iondrag_calc_ions
1043 :
1044 : !=========================================================================
1045 :
1046 0 : subroutine iondrag_calc_ghg (lchnk,ncol,state,ptend)
1047 :
1048 72960 : use phys_grid, only: get_rlat_all_p
1049 : use cam_history, only: outfld
1050 : use physics_types, only: physics_ptend_init
1051 :
1052 :
1053 : !
1054 : ! This subroutine calculates ion drag using globally uniform
1055 : ! ion drag tensor:
1056 : !
1057 : ! |alamxx alamxy |
1058 : ! | |
1059 : ! lambda=| |
1060 : ! | |
1061 : ! |alamyx alamyy |
1062 : !
1063 : ! alamxx and alamxy are provided is data statements
1064 : ! alamyy is obtaine from alamxx:
1065 : !
1066 : !
1067 : ! alamyy = alamxx (sin(DIP_ANGLE))**2
1068 : !
1069 : ! where
1070 : !
1071 : ! DIP_ANGLE = arctan(2.*tan(clat))
1072 : !
1073 :
1074 : !--------------------Input arguments------------------------------------
1075 :
1076 : integer, intent(in) :: lchnk ! chunk identifier
1077 : integer, intent(in) :: ncol ! number of atmospheric columns
1078 :
1079 : type(physics_state), intent(in) :: state
1080 : type(physics_ptend ), intent(out) :: ptend
1081 :
1082 :
1083 :
1084 : !---------------------Local workspace-------------------------------------
1085 :
1086 : real(r8) :: clat(pcols) ! latitudes(radians) for columns
1087 :
1088 : real(r8) alamyyi ! ALAMYY
1089 : real(r8) dipan ! dip angle
1090 : real(r8) dut(pcols,pver)
1091 : real(r8) dvt(pcols,pver)
1092 :
1093 : integer i
1094 : integer k
1095 :
1096 : !-------------------------------------------------------------------------
1097 :
1098 0 : if (.not.doiodrg) then
1099 0 : call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
1100 0 : return
1101 : end if
1102 :
1103 0 : call physics_ptend_init(ptend, state%psetcols, "ion drag", lu=.true., lv=.true.)
1104 :
1105 0 : call get_rlat_all_p(lchnk, pcols, clat)
1106 :
1107 : ! calculate zonal wind drag
1108 0 : dut(:,:)=0.0_r8
1109 0 : do i=1,ncol
1110 0 : do k=ntop_lev,nbot_lev
1111 0 : dut(i,k)=-alamxyi(k)*state%v(i,k)-alamxxi(k)*state%u(i,k)
1112 : enddo
1113 : enddo
1114 :
1115 : ! calculate meridional wind drag
1116 0 : dvt(:,:)=0.0_r8
1117 0 : do i=1,ncol
1118 0 : dipan=atan(2._r8*tan(clat(i)))
1119 0 : do k=ntop_lev,nbot_lev
1120 0 : alamyyi=alamxxi(k)*(sin(dipan))**2._r8
1121 0 : dvt(i,k)=+alamxyi(k)*state%u(i,k)-alamyyi*state%v(i,k)
1122 : enddo
1123 : enddo
1124 :
1125 0 : do i=1,ncol
1126 0 : do k=ntop_lev,nbot_lev
1127 0 : ptend%u(i,k)=dut(i,k)
1128 0 : ptend%v(i,k)=dvt(i,k)
1129 : enddo
1130 : enddo
1131 :
1132 : ! Write out tendencies
1133 0 : call outfld('UIONTEND ',dut ,pcols ,lchnk )
1134 0 : call outfld('VIONTEND ',dvt ,pcols ,lchnk )
1135 :
1136 0 : return
1137 0 : end subroutine iondrag_calc_ghg
1138 :
1139 : !===================================================================================
1140 :
1141 72960 : subroutine iondrag_tend( lchnk, ncol, state, ptend, pbuf, &
1142 : lxx, lyy, lxy, lyx, delt )
1143 :
1144 : !-------------------------------------------------------------------------------
1145 : ! Calculate tendencies in U and V from ion drag tensors, which were
1146 : ! calculated by sub iondrag_calc (module data lxx,lyy,lxy,lyx).
1147 : ! This is called from sub iondrag_calc.
1148 : !-------------------------------------------------------------------------------
1149 :
1150 : !-------------------------------------------------------------------------------
1151 : ! dummy arguments
1152 : !-------------------------------------------------------------------------------
1153 : integer,intent(in) :: lchnk ! current chunk index
1154 : integer,intent(in) :: ncol ! number of atmospheric columns
1155 : real(r8), intent(in) :: delt ! time step (s)
1156 : real(r8), intent(in) :: lxx(pcols,pver) ! ion drag tensor
1157 : real(r8), intent(in) :: lyy(pcols,pver) ! ion drag tensor
1158 : real(r8), intent(in) :: lxy(pcols,pver) ! ion drag tensor
1159 : real(r8), intent(in) :: lyx(pcols,pver) ! ion drag tensor
1160 : type(physics_state), intent(in) :: state ! Physics state variables
1161 : type(physics_ptend), intent(inout) :: ptend ! Physics tendencies
1162 :
1163 : type(physics_buffer_desc), pointer :: pbuf(:)
1164 :
1165 :
1166 : !-------------------------------------------------------------------------------
1167 : ! Local variables
1168 : !-------------------------------------------------------------------------------
1169 : integer :: i, k
1170 : real(r8) :: dti
1171 : real(r8) :: detr
1172 : real(r8) :: us, vs
1173 : real(r8) :: l11, l12, l21, l22
1174 : real(r8) :: dui(pcols,pver) ! zonal ion drag tendency
1175 : real(r8) :: dvi(pcols,pver) ! meridional ion drag tendency
1176 72960 : real(r8), pointer :: ui(:,:) ! pointer to 3d zonal ion drift from edynamo
1177 72960 : real(r8), pointer :: vi(:,:) ! pointer to 3d meridional ion drift from edynamo
1178 72960 : real(r8), pointer :: wi(:,:) ! pointer to 3d vertical ion drift from edynamo
1179 :
1180 : !-------------------------------------------------------------------------------
1181 : ! Get ion ExB drift from physics buffer (they were defined by either the exbdrift
1182 : ! module in chemistry (2d), or the dynamo module in dynamics dpie_coupling (3d),
1183 : ! depending on the switch empirical_ion_velocities. If using dynamo drifts,
1184 : ! they were put into pbuf by dp_coupling. If using empirical exbdrifts, then
1185 : ! they are redundant in the vertical dimension (i.e., 2d only).
1186 : !-------------------------------------------------------------------------------
1187 72960 : call pbuf_get_field(pbuf, ui_idx, ui)
1188 72960 : call pbuf_get_field(pbuf, vi_idx, vi)
1189 72960 : call pbuf_get_field(pbuf, wi_idx, wi)
1190 :
1191 72960 : dti = 1._r8/delt
1192 : !-------------------------------------------------------------------------------
1193 : ! Zonal (du) and meridional (dv) wind drag, using ExB drift velocities
1194 : ! from exbdrift module (pbuf):
1195 : !-------------------------------------------------------------------------------
1196 1824000 : do k = ntop_lev,nbot_lev
1197 27038976 : do i = 1,ncol
1198 : !-------------------------------------------------------------------------------
1199 : ! 2/28/04 btf:
1200 : ! Full ion-drag, using lambdas and ExB drifts.
1201 : ! This should succeed with bz = 0 (efield module)
1202 : ! Runs:
1203 : ! bz=-5, nstep=24 min, nsplit=4 (6 min dynamics): crashed after 2 days.
1204 : ! bz=-5, nstep=24 min, nsplit=6 (4 min dynamics): 5 day run succeeded.
1205 : ! See comments in efield module re bz < 0 (efield.F90).
1206 : ! Ion drifts are from edynamo if use_dynamo_drifts=true, exbdrift otherwise.
1207 : !-------------------------------------------------------------------------------
1208 25214976 : us = ui(i,k) - state%u(i,k)
1209 25214976 : vs = vi(i,k) - state%v(i,k)
1210 :
1211 : !-------------------------------------------------------------------------------
1212 : ! Exclude ue,ve drift momentum source to avoid crashes when bz < 0 and
1213 : ! full 30 min timestep (partial ion-drag):
1214 : !-------------------------------------------------------------------------------
1215 25214976 : l11 = dti + lxx(i,k)
1216 25214976 : l12 = lxy(i,k)
1217 25214976 : l21 = -lyx(i,k)
1218 25214976 : l22 = dti + lyy(i,k)
1219 25214976 : detr = dti/(l11*l22 - l12*l21)
1220 25214976 : dui(i,k) = dti*(detr*(l12*vs - l22*us) + us)
1221 26966016 : dvi(i,k) = dti*(detr*(l21*us - l11*vs) + vs)
1222 : end do
1223 : end do
1224 :
1225 : !-------------------------------------------------------------------------------
1226 : ! Apply to model tendencies:
1227 : !-------------------------------------------------------------------------------
1228 1824000 : do k = ntop_lev,nbot_lev
1229 : !-------------------------------------------------------------------------------
1230 : ! Ion drag tendencies:
1231 : !-------------------------------------------------------------------------------
1232 26966016 : ptend%u(:ncol,k) = dui(:ncol,k)
1233 27038976 : ptend%v(:ncol,k) = dvi(:ncol,k)
1234 : !-------------------------------------------------------------------------------
1235 : ! Turn off ion drag tendency:
1236 : !-------------------------------------------------------------------------------
1237 : ! ptend%u(:ncol,k) = 0._r8
1238 : ! ptend%v(:ncol,k) = 0._r8
1239 : end do
1240 3429120 : do k = nbot_lev+1,pver
1241 51684864 : dui(:ncol,k) = 0._r8
1242 51684864 : dvi(:ncol,k) = 0._r8
1243 51684864 : ptend%u(:ncol,k) = 0._r8
1244 51757824 : ptend%v(:ncol,k) = 0._r8
1245 : end do
1246 :
1247 : !
1248 : ! Ion drifts are either empirical (use_dynamo_drifts==false), or from edynamo
1249 : ! (use_dynamo_drifts==true). See addfld calls in this source file.
1250 : ! If empirical, the drifts will be 2d (i.e., redundant in the vertical dimension)
1251 : !
1252 72960 : if (empirical_ion_velocities) then
1253 72960 : call outfld ( 'UI', ui, pcols, lchnk )
1254 72960 : call outfld ( 'VI', vi, pcols, lchnk )
1255 72960 : call outfld ( 'WI', wi, pcols, lchnk )
1256 : endif
1257 :
1258 72960 : call outfld ( 'UIONTEND', dui, pcols, lchnk ) ! u ion drag tendency
1259 72960 : call outfld ( 'VIONTEND', dvi, pcols, lchnk ) ! v ion drag tendency
1260 :
1261 72960 : end subroutine iondrag_tend
1262 :
1263 : !================================================================================================
1264 72960 : subroutine jouleheat_tend( lchnk, ncol, state, ptend, pbuf, &
1265 : lxx, lyy, lxy, lyx )
1266 : !-------------------------------------------------------------------------------
1267 : ! Calculate tendencies in T due to joule heating.
1268 : ! This is called from sub iondrag_calc.
1269 : !-------------------------------------------------------------------------------
1270 :
1271 : use physconst, only: pi
1272 : use air_composition, only: cpairv
1273 : use phys_grid, only: get_rlon_p, get_rlat_p
1274 :
1275 : !-------------------------------------------------------------------------------
1276 : ! dummy arguments
1277 : !-------------------------------------------------------------------------------
1278 : integer,intent(in) :: lchnk ! current chunk index
1279 : integer,intent(in) :: ncol ! number of atmospheric columns
1280 : real(r8), intent(in) :: lxx(pcols,pver) ! ion drag tensor
1281 : real(r8), intent(in) :: lyy(pcols,pver) ! ion drag tensor
1282 : real(r8), intent(in) :: lxy(pcols,pver) ! ion drag tensor
1283 : real(r8), intent(in) :: lyx(pcols,pver) ! ion drag tensor
1284 : type(physics_state), intent(in) :: state ! Physics state variables
1285 : type(physics_ptend), intent(inout) :: ptend ! Physics tendencies (inout)
1286 :
1287 : type(physics_buffer_desc), pointer :: pbuf(:)
1288 :
1289 : !-------------------------------------------------------------------------------
1290 : ! Local variables
1291 : !-------------------------------------------------------------------------------
1292 : integer :: k, i
1293 : integer :: max_ind(2)
1294 : real(r8) :: us, vs
1295 : real(r8) :: max_q
1296 : real(r8) :: qjoule(pcols,pver) ! joule heating
1297 : real(r8) :: qout(pcols,pver) ! temp for outfld
1298 72960 : real(r8), pointer :: ui(:,:) ! pointer to pbuf
1299 72960 : real(r8), pointer :: vi(:,:) ! pointer to pbuf
1300 :
1301 : logical, parameter :: debug = .false.
1302 :
1303 : !-------------------------------------------------------------------------------
1304 : ! Get ion velocities from physics buffer (they were defined by exbdrift module)
1305 : ! Ion velocities are 2d arrays, i.e., no vertical dimension.
1306 : !-------------------------------------------------------------------------------
1307 72960 : call pbuf_get_field(pbuf, ui_idx, ui )
1308 72960 : call pbuf_get_field(pbuf, vi_idx, vi )
1309 :
1310 1824000 : do k = ntop_lev,nbot_lev
1311 : ! write(iulog,"('qjoule: k=',i3,' u=',/,(6e12.4))") k,state%u(:,k)
1312 : ! write(iulog,"('qjoule: k=',i3,' v=',/,(6e12.4))") k,state%v(:,k)
1313 : ! write(iulog,"('qjoule: k=',i3,' lxx=',/,(6e12.4))") k,lxx(:,k)
1314 : ! write(iulog,"('qjoule: k=',i3,' lxy=',/,(6e12.4))") k,lxy(:,k)
1315 : ! write(iulog,"('qjoule: k=',i3,' lyx=',/,(6e12.4))") k,lyx(:,k)
1316 : ! write(iulog,"('qjoule: k=',i3,' lyy=',/,(6e12.4))") k,lyy(:,k)
1317 27038976 : do i = 1,ncol
1318 25214976 : us = ui(i,k) - state%u(i,k)
1319 25214976 : vs = vi(i,k) - state%v(i,k)
1320 25214976 : qjoule(i,k) = us*us*lxx(i,k) + us*vs*(lxy(i,k) - lyx(i,k)) + vs*vs*lyy(i,k)
1321 26966016 : ptend%s(i,k) = qjoule(i,k) ! joule heating tendency
1322 : end do
1323 : ! write(iulog,"('qjoule: k=',i3,' qjoule(:,k)=',/,(6e12.4))") k,qjoule(:,k)
1324 : end do
1325 3429120 : do k = nbot_lev+1,pver
1326 51684864 : qjoule(:ncol,k) = 0._r8
1327 51757824 : ptend%s(:ncol,k) = 0._r8 ! no joule heating tendency
1328 : end do
1329 :
1330 : sw_debug: if (debug) then
1331 : max_q = 100._r8*maxval( abs( qjoule(:ncol,ntop_lev:nbot_lev) )/state%t(:ncol,ntop_lev:nbot_lev) )
1332 : max_ind(:) = maxloc( abs( qjoule(:ncol,ntop_lev:nbot_lev) )/state%t(:ncol,ntop_lev:nbot_lev) )
1333 : i = max_ind(1)
1334 : k = max_ind(2)
1335 : if( lchnk == 25 ) then
1336 : i = 14
1337 : k = 3
1338 : write(iulog,*) ' '
1339 : write(iulog,*) '-------------------------------------------------------'
1340 : write(iulog,*) 'jouleheat_tend: lon,lat = ',get_rlon_p(lchnk,14)*180._r8/pi, get_rlat_p(lchnk,14)*180._r8/pi
1341 : write(iulog,*) 'jouleheat_tend: dt,t,max% dt/t = ',qjoule(i,k)/cpairv(i,k,lchnk),state%t(i,k),max_q, &
1342 : ' @ lchnk,i,k = ',lchnk,max_ind(:)
1343 : write(iulog,*) 'jouleheat_tend: lxx,xy,yx,yy = ',lxx(i,k),lxy(i,k),lyx(i,k),lyy(i,k)
1344 : write(iulog,*) 'jouleheat_tend: u,ui,v,vi = ',state%u(i,k),ui(i,k),state%v(i,k),vi(i,k)
1345 : write(iulog,*) 'jouleheat_tend: us,vs = ',ui(i,k) - state%u(i,k),vi(i,k) - state%v(i,k)
1346 : write(iulog,*) 'jouleheat_tend: du,dv = ',ptend%u(i,k),ptend%v(i,k)
1347 : write(iulog,*) 'jouleheat_tend: dt'
1348 : write(iulog,'(1p,5g15.7)') qjoule(max_ind(1),ntop_lev:nbot_lev)/cpairv(max_ind(1),ntop_lev:nbot_lev,lchnk)
1349 : write(iulog,*) '-------------------------------------------------------'
1350 : write(iulog,*) ' '
1351 : ! stop 'diagnostics'
1352 : end if
1353 : endif sw_debug
1354 :
1355 78723840 : qout(:ncol,:) = qjoule(:ncol,:)/cpairv(:ncol,:,lchnk)
1356 :
1357 72960 : call outfld ( 'QJOULE', qout, pcols, lchnk )
1358 :
1359 145920 : end subroutine jouleheat_tend
1360 :
1361 : !==============================================================================
1362 :
1363 0 : subroutine iondrag_inidat(ncid_ini, pbuf2d)
1364 :
1365 72960 : use pio, only: file_desc_t
1366 : use ncdio_atm,only: infld
1367 : use infnan, only: nan, assignment(=)
1368 : use cam_grid_support, only : cam_grid_check, cam_grid_id, cam_grid_get_dim_names
1369 : use physics_buffer, only : pbuf_set_field
1370 :
1371 : ! args
1372 : type(file_desc_t), intent(inout) :: ncid_ini ! Initial condition file id
1373 : type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer
1374 :
1375 : ! local vars
1376 0 : real(r8), pointer :: ui_tmp(:,:,:)
1377 0 : real(r8), pointer :: vi_tmp(:,:,:)
1378 0 : real(r8), pointer :: wi_tmp(:,:,:)
1379 : real(r8) :: nanval
1380 : integer :: grid_id
1381 : character(len=4) :: dim1name, dim2name
1382 : character(len=*), parameter :: subname='iondrag_inidat'
1383 : logical :: found_ui, found_vi, found_wi
1384 :
1385 0 : allocate(ui_tmp(pcols,pver,begchunk:endchunk))
1386 0 : allocate(vi_tmp(pcols,pver,begchunk:endchunk))
1387 0 : allocate(wi_tmp(pcols,pver,begchunk:endchunk))
1388 :
1389 0 : grid_id = cam_grid_id('physgrid')
1390 0 : if (.not. cam_grid_check(grid_id)) then
1391 0 : call endrun(trim(subname)//': Internal error, no "physgrid" grid')
1392 : end if
1393 0 : call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
1394 :
1395 : call infld( 'UI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
1396 0 : ui_tmp, found_ui, gridname='physgrid')
1397 : call infld( 'VI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
1398 0 : vi_tmp, found_vi, gridname='physgrid')
1399 : call infld( 'WI',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
1400 0 : wi_tmp, found_wi, gridname='physgrid')
1401 :
1402 0 : ionvels_read_from_file = found_ui .and. found_vi .and. found_wi
1403 :
1404 0 : ui_idx = pbuf_get_index('UI')
1405 0 : vi_idx = pbuf_get_index('VI')
1406 0 : wi_idx = pbuf_get_index('WI')
1407 :
1408 0 : if (ionvels_read_from_file) then
1409 0 : call pbuf_set_field(pbuf2d, ui_idx, ui_tmp)
1410 0 : call pbuf_set_field(pbuf2d, vi_idx, vi_tmp)
1411 0 : call pbuf_set_field(pbuf2d, wi_idx, wi_tmp)
1412 : else
1413 0 : nanval=nan
1414 0 : call pbuf_set_field(pbuf2d, ui_idx, nanval)
1415 0 : call pbuf_set_field(pbuf2d, vi_idx, nanval)
1416 0 : call pbuf_set_field(pbuf2d, wi_idx, nanval)
1417 : endif
1418 :
1419 0 : deallocate( ui_tmp )
1420 0 : deallocate( vi_tmp )
1421 0 : deallocate( wi_tmp )
1422 :
1423 0 : if (masterproc) then
1424 0 : write(iulog,*) 'iondrag_inidat: ionvels_read_from_file = ',ionvels_read_from_file
1425 : end if
1426 :
1427 0 : end subroutine iondrag_inidat
1428 :
1429 : end module iondrag
|
