Line data Source code
1 : module ion_electron_temp
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : !
6 : ! Module to compute the ion/electron temperature and dry static heating
7 : !
8 : ! Authors: Joe McInerney/Hanli Liu/Art Richmond
9 : !
10 : !---------------------------------------------------------------------------------
11 : use shr_kind_mod, only : r8 => shr_kind_r8 ! Real kind to declare variables
12 : use ppgrid, only : pcols, pver, pverp ! Dimensions and chunk bounds
13 : use ppgrid, only : begchunk, endchunk
14 : use cam_history, only : outfld, hist_fld_active, write_inithist ! Routine to output fields to history files
15 : use cam_history, only : horiz_only, addfld, add_default ! Routines and variables for adding fields to history output
16 : use physics_types, only : physics_state, & ! Structures containing physics state variables
17 : physics_ptend, & ! Structures containing physics tendency variables
18 : physics_ptend_init ! Routine to initialize physics tendency variables
19 : use physics_buffer, only : pbuf_add_field, pbuf_get_chunk, &
20 : pbuf_get_index,dtype_r8, & !
21 : physics_buffer_desc, & !
22 : pbuf_get_field, & ! Needed to access physics buffer
23 : pbuf_set_field
24 : use mo_jeuv, only : nIonRates ! Number of ionization rates in mo_photo
25 : use shr_const_mod, only : kboltz => shr_const_boltz, &
26 : pi => shr_const_pi ! Boltzmann constant and pi
27 : use chem_mods, only : adv_mass ! Array holding mass values for short lived species
28 : use cam_abortutils, only : endrun
29 : use mo_chem_utls, only : get_spc_ndx ! Routine to get index of adv_mass array for short lived species
30 : use constituents, only : cnst_get_ind, cnst_mw ! Routines to get molecular weights for constituents
31 : use solar_parms_data, only : f107=>solar_parms_f107 ! 10.7 cm solar flux
32 : use steady_state_tei, only : steady_state_tei_init, steady_state_tei_tend
33 : use perf_mod, only : t_startf, t_stopf ! timing utils
34 : use spmd_utils, only : masterproc
35 : use cam_logfile, only : iulog ! Output unit
36 : use ionos_state_mod, only : ionos_state
37 : use air_composition,only : cpairv
38 :
39 : implicit none
40 :
41 : save
42 :
43 : private ! Make default type private to the module
44 :
45 : !------------------------
46 : ! PUBLIC: interfaces
47 : !------------------------
48 : public :: ion_electron_temp_init ! Initialization
49 : public :: ion_electron_temp_timestep_init
50 : public :: ion_electron_temp_register ! Registration of ionosphere variables in pbuf physics buffer
51 : public :: ion_electron_temp_inidat ! Get fields from initial condition file into physics buffer
52 : public :: ion_electron_temp_tend ! Calculate tendencies for extended model ionosphere
53 : public :: ion_electron_temp_readnl
54 :
55 : !------------------------------------------------------------------------
56 : ! PRIVATE: Rest of the data and interfaces are private to this module
57 : !------------------------------------------------------------------------
58 : real(r8), parameter :: kboltz_ev = 8.617E-5_r8 ! Boltzmann constant (eV/K)
59 : real(r8), parameter :: temax = 7.0E3_r8 ! maximum electron temperature (K)
60 : real(r8), parameter :: dayOPFlux = 2.0E8_r8 ! Daytime O+ flux at upper boundary (
61 : real(r8), parameter :: nightOPFlux = -2.0E8_r8 ! Nighttime O+ flux at upper boundary (
62 :
63 : real(r8), parameter :: rads2Degs = 180._r8/pi ! radians to degrees
64 :
65 : ! private data
66 : real(r8) :: rMassOp ! O+ molecular weight kg/kmol
67 :
68 : logical :: steady_state_ion_elec_temp = .true.
69 : logical :: initialized_TiTe = .false.
70 :
71 : integer :: index_te=-1, index_ti=-1 ! Indices to find ion and electron temperature in pbuf
72 :
73 : contains
74 :
75 : !==============================================================================
76 :
77 768 : subroutine ion_electron_temp_readnl(nlfile)
78 :
79 : use namelist_utils, only: find_group_name
80 : use units, only: getunit, freeunit
81 : use spmd_utils, only: mpicom, masterprocid, mpicom, mpi_logical
82 :
83 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
84 :
85 : ! Local variables
86 : integer :: unitn, ierr
87 : character(len=*), parameter :: subname = 'ion_electron_temp_readnl'
88 :
89 : namelist /ion_electron_temp_nl/ steady_state_ion_elec_temp
90 :
91 768 : if (masterproc) then
92 2 : unitn = getunit()
93 2 : open( unitn, file=trim(nlfile), status='old' )
94 2 : call find_group_name(unitn, 'ion_electron_temp_nl', status=ierr)
95 2 : if (ierr == 0) then
96 2 : read(unitn, ion_electron_temp_nl, iostat=ierr)
97 2 : if (ierr /= 0) then
98 0 : call endrun(subname // ':: ERROR reading namelist')
99 : end if
100 : end if
101 2 : close(unitn)
102 2 : call freeunit(unitn)
103 :
104 : end if
105 :
106 768 : call mpi_bcast (steady_state_ion_elec_temp, 1, mpi_logical, masterprocid, mpicom, ierr)
107 :
108 768 : if (masterproc) then
109 2 : write(iulog,*) subname//': steady_state_ion_elec_temp = ',steady_state_ion_elec_temp
110 : endif
111 :
112 768 : end subroutine ion_electron_temp_readnl
113 :
114 : !==============================================================================
115 : !==============================================================================
116 :
117 1536 : subroutine ion_electron_temp_init(pbuf2d)
118 :
119 : !-----------------------------------------------------------------------
120 : ! Time independent initialization for ionosphere simulation.
121 : !-----------------------------------------------------------------------
122 :
123 : use phys_control, only: phys_getopts !Method used to get flag for waccmx ionosphere output variables
124 :
125 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
126 :
127 : logical :: history_waccmx
128 : integer :: indxOp,sIndxOp ! state%q or pbuf index for O+ mixing ratio
129 :
130 768 : if (steady_state_ion_elec_temp) then
131 0 : call steady_state_tei_init(pbuf2d)
132 : end if
133 :
134 768 : call phys_getopts(history_waccmx_out=history_waccmx)
135 :
136 : !-------------------------------------------------------------------------------
137 : ! Add history variables for ionosphere
138 : !-------------------------------------------------------------------------------
139 1536 : call addfld ('QIonElec' ,(/ 'lev' /), 'I', 'K sec-1', 'Electron Ion Thermal Heating Rate')
140 1536 : call addfld ('TElec&IC' ,(/ 'lev' /), 'I', 'K', 'Electron Temperature')
141 1536 : call addfld ('TIon&IC' ,(/ 'lev' /), 'I', 'K', 'Ion Temperature')
142 1536 : call addfld ('TElec' ,(/ 'lev' /), 'I', 'K', 'Electron Temperature')
143 1536 : call addfld ('TIon' ,(/ 'lev' /), 'I', 'K', 'Ion Temperature')
144 768 : call addfld ('ElecColDens' ,horiz_only , 'I', 'TECU', 'Electron Column Density')
145 768 : if (.not.steady_state_ion_elec_temp) then
146 1536 : call addfld ('QIN' ,(/ 'lev' /), 'I', 'K sec-1','Ion-neutral Heating Rate')
147 1536 : call addfld ('QEN' ,(/ 'lev' /), 'I', 'K sec-1','Electron-neutral Heating Rate')
148 1536 : call addfld ('QEI' ,(/ 'lev' /), 'I', 'K sec-1','Electron-ion Heating Rate')
149 1536 : call addfld ('LOSS_g3' ,(/ 'lev' /), 'I', ' ', 'Loss Term g3')
150 1536 : call addfld ('LOSS_EI' ,(/ 'lev' /), 'I', ' ', 'Loss Term EI')
151 1536 : call addfld ('LOSS_IN' ,(/ 'lev' /), 'I', ' ', 'Loss Term IN')
152 1536 : call addfld ('SOURCER' ,(/ 'lev' /), 'I', ' ', 'SOURCER')
153 1536 : call addfld ('SOURCEEff' ,(/ 'lev' /), 'I', ' ', 'SOURCEEff')
154 1536 : call addfld ('AURIPRATESUM' ,(/ 'lev' /), 'I', ' ', 'Auroral ionization')
155 1536 : call addfld ('OpI' ,(/ 'lev' /), 'I', ' ', 'O+ Ionosphere')
156 1536 : call addfld ('eI' ,(/ 'lev' /), 'I', ' ', 'e Ionosphere')
157 : end if
158 :
159 768 : call add_default ('TElec&IC' , 0, ' ')
160 768 : call add_default ('TIon&IC' , 0, ' ')
161 :
162 : !-------------------------------------------------------------------------------
163 : ! Set default values for ionosphere history variables
164 : !-------------------------------------------------------------------------------
165 768 : if (history_waccmx) then
166 768 : call add_default ('TElec' , 1, ' ')
167 768 : call add_default ('TIon' , 1, ' ')
168 768 : if (.not.steady_state_ion_elec_temp) then
169 768 : call add_default ('QIN' , 1, ' ')
170 768 : call add_default ('QEN' , 1, ' ')
171 768 : call add_default ('QEI' , 1, ' ')
172 768 : call add_default ('SOURCER' , 1, ' ')
173 768 : call add_default ('SOURCEEff' , 1, ' ')
174 768 : call add_default ('AURIPRATESUM' , 1, ' ')
175 : end if
176 : end if
177 :
178 768 : if (.not.steady_state_ion_elec_temp) then
179 768 : call cnst_get_ind( 'Op', indxOp, abort=.false. )
180 768 : if (indxOp > 0) then
181 0 : rMassOp = cnst_mw(indxOP)
182 : else
183 768 : sIndxOp = get_spc_ndx( 'Op' )
184 768 : if (sIndxOp > 0) then
185 768 : rMassOp = adv_mass(sIndxOp)
186 : else
187 0 : call endrun('update_teti: Cannot find short-lived index for Op in update_teti')
188 : endif
189 : endif
190 : endif
191 :
192 768 : end subroutine ion_electron_temp_init
193 :
194 : !==============================================================================
195 16128 : subroutine ion_electron_temp_timestep_init(phys_state,pbuf2d)
196 : use time_manager, only: is_first_step
197 :
198 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
199 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
200 :
201 : integer :: ncol
202 : integer :: lchnk
203 8064 : type(physics_buffer_desc), pointer :: phys_buffer_chunk(:)
204 :
205 8064 : if (is_first_step() .and. .not.initialized_TiTe .and. index_te>0 .and. index_ti>0) then
206 0 : do lchnk=begchunk,endchunk
207 0 : ncol=phys_state(lchnk)%ncol
208 0 : phys_buffer_chunk => pbuf_get_chunk(pbuf2d,lchnk)
209 0 : call pbuf_set_field(phys_buffer_chunk,index_te,phys_state(lchnk)%t(:ncol,:),start=(/1,1/), kount=(/ncol,pver/))
210 0 : call pbuf_set_field(phys_buffer_chunk,index_ti,phys_state(lchnk)%t(:ncol,:),start=(/1,1/), kount=(/ncol,pver/))
211 : end do
212 0 : initialized_TiTe=.true.
213 : endif
214 16128 : end subroutine ion_electron_temp_timestep_init
215 :
216 : !==============================================================================
217 :
218 768 : subroutine ion_electron_temp_register
219 :
220 : !-----------------------------------------------------------------------
221 : ! Register ionosphere variables with physics buffer:
222 : !
223 : ! Ion production rates pcols,pver,nIonRates,
224 : ! so firstdim = 1 middledim = pver lastdim = nIonRates.
225 : !
226 : ! pcols dimension and lchnk assumed here
227 : !
228 : !-----------------------------------------------------------------------
229 :
230 : !------------------------------------------------------------------------------
231 : ! Electron temperature in physics buffer (global so can write to history files)
232 : !------------------------------------------------------------------------------
233 768 : call pbuf_add_field('TElec','global',dtype_r8,(/pcols,pver/), index_te)
234 :
235 : !--------------------------------------------------------------------------
236 : ! Ion temperature in physics buffer (global so can write to history files)
237 : !--------------------------------------------------------------------------
238 768 : call pbuf_add_field('TIon', 'global',dtype_r8,(/pcols,pver/), index_ti)
239 :
240 8064 : end subroutine ion_electron_temp_register
241 :
242 : !==============================================================================
243 :
244 384 : subroutine ion_electron_temp_inidat(ncid_ini, pbuf2d)
245 :
246 : !-----------------------------------------------------------------------
247 : ! Grab fields from initial condition file and put in physics buffer
248 : !-----------------------------------------------------------------------
249 :
250 : use pio, only: file_desc_t
251 : use cam_grid_support, only: cam_grid_check, cam_grid_id
252 : use cam_grid_support, only: cam_grid_get_dim_names
253 : use cam_abortutils, only: endrun
254 : use ncdio_atm, only: infld
255 : use ppgrid, only: pcols, pver, begchunk, endchunk
256 : use infnan, only: nan, assignment(=)
257 :
258 : type(file_desc_t), intent(inout) :: ncid_ini ! Initial condition file id
259 : type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer
260 :
261 : integer :: grid_id
262 : character(len=4) :: dim1name, dim2name
263 : logical :: found
264 384 : real(r8),pointer :: tE(:,:,:) ! Electron temperature pointer
265 384 : real(r8),pointer :: tI(:,:,:) ! Ion temperature pointer
266 : integer :: ierr
267 : character(len=*), parameter :: subname='ION_ELECTRON_TEMP_INIDAT'
268 : real(r8) :: nanval
269 :
270 384 : nanval=nan
271 384 : found = .false.
272 :
273 384 : grid_id = cam_grid_id('physgrid')
274 384 : if (.not. cam_grid_check(grid_id)) then
275 0 : call endrun(trim(subname)//': Internal error, no "physgrid" grid')
276 : end if
277 384 : call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
278 :
279 384 : if (index_te>0) then
280 : !---------------------------------------------------------------------------------
281 : ! Electron temperature
282 : !---------------------------------------------------------------------------------
283 1152 : allocate(tE(pcols,pver,begchunk:endchunk))
284 : call infld( 'TElec',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
285 384 : tE, found, gridname='physgrid')
286 :
287 384 : if (.not.found) then
288 0 : if (masterproc) write(iulog,*) 'ion_electron_temp_inidat: Could not find electron temperature in ic file. ' &
289 0 : // 'Try to read neutral temperature.'
290 : call infld( 'T',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
291 0 : tE, found, gridname='physgrid')
292 : endif
293 384 : if (found) then
294 384 : call pbuf_set_field(pbuf2d, index_te, tE)
295 : else
296 0 : call pbuf_set_field(pbuf2d, index_te, nanval)
297 : endif
298 :
299 384 : initialized_TiTe = found
300 384 : deallocate(tE)
301 : endif
302 :
303 384 : if (index_ti>0) then
304 : !----------------------------------------------------------------------------
305 : ! Ion temperature
306 : !----------------------------------------------------------------------------
307 1152 : allocate(tI(pcols,pver,begchunk:endchunk))
308 384 : if (initialized_TiTe) then ! try to initialize ion temp only if electron temp was initialized above
309 : call infld( 'TIon',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
310 384 : tI, found, gridname='physgrid')
311 384 : if (.not.found) then
312 0 : if (masterproc) write(iulog,*) 'ion_electron_temp_inidat: Could not find ion temperature in ic file. ' &
313 0 : // 'Try to read neutral temperature.'
314 : call infld( 'T',ncid_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, &
315 0 : tI, found, gridname='physgrid')
316 : endif
317 : endif
318 384 : if (found) then
319 384 : call pbuf_set_field(pbuf2d, index_ti, tI)
320 : else
321 0 : call pbuf_set_field(pbuf2d, index_ti, nanval)
322 : endif
323 :
324 384 : initialized_TiTe = initialized_TiTe .and. found
325 384 : deallocate(tI)
326 : endif
327 384 : if (index_te>0 .and. index_ti>0 .and. .not.initialized_TiTe) then
328 : write(iulog,*) 'ion_electron_temp_inidat: Not able to read temperatures from IC file.' &
329 0 : //' Will set ion and electron temperatures to neutral temperature (state%t) on initial timestep.'
330 : endif
331 :
332 768 : end subroutine ion_electron_temp_inidat
333 :
334 : !==============================================================================
335 :
336 21888 : subroutine ion_electron_temp_tend(state, ptend, pbuf, ztodt)
337 :
338 : !-------------------------------------------------------------------------------------
339 : ! Calculate dry static energy and O+ tendency for extended ionosphere simulation
340 : !-------------------------------------------------------------------------------------
341 :
342 : !------------------------------Arguments--------------------------------
343 :
344 384 : use physics_types, only: physics_ptend_sum
345 :
346 : type(physics_state), intent(in) :: state ! physics state structure
347 : type(physics_ptend), intent(inout) :: ptend ! parameterization tendency structure
348 : type(physics_buffer_desc),pointer :: pbuf(:) ! physics buffer
349 :
350 : real(r8), intent(in) :: ztodt ! Physics time step
351 :
352 : !---------------------------Local storage-------------------------------
353 :
354 2232576 : type(physics_ptend) :: ptend_loc ! Local parameterization tendencies
355 : type(ionos_state) :: istate ! ionosphere state structure
356 :
357 : integer :: lchnk ! Chunk number
358 : integer :: ncol ! Number of columns in chunk
359 :
360 : integer :: teTiBot ! bottom of ionosphere calculations
361 :
362 21888 : real(r8), dimension(:,:), pointer :: tE ! Pointer to electron temperature in pbuf (K)
363 21888 : real(r8), dimension(:,:), pointer :: tI ! Pointer to ion temperature in pbuf (K)
364 :
365 : logical :: ls
366 : real(r8) :: dse_tend(pcols,pver) ! dry static energy tendency
367 : real(r8) :: qionelec(pcols,pver) ! diagnostic heating rate (neutrals)
368 :
369 21888 : call t_startf ('ion_electron_temp_tend')
370 :
371 : !----------------------------------------------------------------
372 : ! Get number of this chunk
373 : !----------------------------------------------------------------
374 21888 : lchnk = state%lchnk
375 21888 : ncol = state%ncol
376 :
377 21888 : ls = .TRUE.
378 21888 : call physics_ptend_init(ptend_loc, state%psetcols, 'ionosphere', ls=ls)
379 :
380 : !-------------------------------------------------------------------------------------------------------------------
381 : ! Get electron and ion temperatures from physics buffer.
382 : !-------------------------------------------------------------------------------------------------------------------
383 21888 : call pbuf_get_field(pbuf, index_te, tE)
384 21888 : call pbuf_get_field(pbuf, index_ti, tI)
385 :
386 : !------------------------------------------------------------
387 : ! Initialize data needed in the ionosphere calculations
388 : !------------------------------------------------------------
389 21888 : call update_istate(state, pbuf, istate, teTiBot)
390 :
391 21888 : if (steady_state_ion_elec_temp) then
392 : !-----------------------------------------------------------------
393 : ! steady-state solution
394 : !-----------------------------------------------------------------
395 0 : dse_tend(:ncol,:) = ptend%s(:ncol,:)
396 0 : call steady_state_tei_tend(state, istate, dse_tend, pbuf)
397 0 : ptend_loc%s(:ncol,:) = dse_tend(:ncol,:)
398 : else
399 : !-----------------------------------------------------------------
400 : ! Get electron temperature and update dry static energy tendency
401 : !-----------------------------------------------------------------
402 21888 : call update_teti(state, ptend%s, ptend_loc%s, ztodt, istate, tE, tI, teTiBot)
403 : end if
404 :
405 : !--------------------------------------------------------------
406 : ! Make Te and Ti fields available for output to history files
407 : !--------------------------------------------------------------
408 21888 : call outfld ('TElec' , tE, pcols, lchnk)
409 21888 : call outfld ('TIon' , tI, pcols, lchnk)
410 21888 : if (write_inithist()) then
411 0 : call outfld ('TElec&IC', tE, pcols, lchnk)
412 0 : call outfld ('TIon&IC' , tI, pcols, lchnk)
413 : endif
414 :
415 37012608 : qionelec(:ncol,:) = ptend_loc%s(:ncol,:)/cpairv(:ncol,:,lchnk)
416 21888 : call outfld ('QIonElec' , qionelec, pcols, lchnk)
417 :
418 21888 : call physics_ptend_sum(ptend_loc, ptend, ncol)
419 :
420 21888 : call t_stopf ('ion_electron_temp_tend')
421 :
422 43776 : end subroutine ion_electron_temp_tend
423 :
424 : !===============================================================================
425 :
426 21888 : subroutine update_istate(state, pbuf, istate, teTiBot)
427 :
428 : !---------------------------------------------------------------------------------------
429 : ! Time independent initialization for extended ionosphere simulation called in phys_init
430 : ! of physpkg module which is called in cam_comp module
431 : !---------------------------------------------------------------------------------------
432 21888 : use mo_apex, only: bnorth, beast, bdown ! Magnetic field components
433 : use time_manager, only: get_curr_calday ! Routine to get current calendar day
434 : use physconst, only: rearth
435 : use air_composition, only: rairv, mbarv ! Constituent dependent rair and mbar
436 : use ref_pres, only: press_lim_idx
437 : use orbit, only: zenith
438 :
439 : use short_lived_species, only: slvd_index,slvd_pbf_ndx => pbf_idx ! Routines to access short lived species
440 :
441 : type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer
442 : type(physics_state), intent(in), target :: state ! physics state structure
443 : type(ionos_state), intent(inout), target :: istate ! ionosphere state structure
444 :
445 : integer, intent(out) :: teTiBot ! bottom of ionosphere calculations
446 :
447 : !---------------------------Local storage-------------------------------
448 : integer,parameter :: nCnst = 9 ! Number of species needed from state%q or pbuf
449 :
450 : integer :: lchnk ! Chunk number
451 : integer :: ncol ! Number of columns in current chunk
452 :
453 : integer :: indxIR ! pbuf index for ionization rates
454 : integer :: indxAIPRS ! pbuf index for aurora ion production rate sum
455 :
456 : integer :: indxCnst ! Constituent index used in cslculating densities
457 :
458 : integer :: indxSLvd ! index of pbuf to access short lived species
459 : integer :: sIndx ! index of adv_mass for any short lived species to access constituent mass
460 :
461 : integer :: iVer ! Counter for vertical loops
462 : integer :: iCol ! Counter for column loops
463 : integer :: iIonR ! Counter for ionization rates loops
464 : integer :: iCnst ! Counter for constituent loop
465 :
466 : integer :: indxSP ! pbuf index for Pedersen Conductivity
467 : integer :: indxSH ! pbuf index for Hall Conductivity
468 :
469 : real(r8), parameter :: teTiBotPres = 50._r8 ! Pressure above which electron/ion temperature are calculated
470 : ! in WACCM-X. (Pa)
471 :
472 : character(len = 3), dimension(nCnst) :: cCnst
473 :
474 21888 : real(r8), dimension(:,:), pointer :: sigma_ped ! Pointer to Pedersen Conductivity in pbuf (siemens/m) from module iondrag
475 21888 : real(r8), dimension(:,:), pointer :: sigma_hall ! Pointer to Hall Conductivity in pbuf (siemens/m)
476 :
477 21888 : real(r8), dimension(:,:), pointer :: mmrP ! Pointer to access short lived species in pbuf
478 :
479 21888 : real(r8), dimension(:),pointer :: geoLatR ! Latitude (radians) Make ncol because zenith aurora are ncol
480 21888 : real(r8), dimension(:),pointer :: geoLonR ! Longitude (radians)
481 :
482 21888 : real(r8), dimension(:,:),pointer :: pMid ! Midpoint pressure (Pa)
483 21888 : real(r8), dimension(:,:),pointer :: tN ! Neutral temperature (K)
484 :
485 21888 : real(r8), dimension(:,:),pointer :: tNInt ! Interface Temperture (K)
486 :
487 21888 : real(r8), dimension(:),pointer :: cosZenAngR ! cosine of zenith angle (radians)
488 21888 : real(r8), dimension(:),pointer :: zenAngD ! zenith angle (degrees)
489 :
490 21888 : real(r8), dimension(:,:),pointer :: bNorth3d ! northward component of magnetic field units?
491 21888 : real(r8), dimension(:,:),pointer :: bEast3d ! eastward component of magnetic field
492 21888 : real(r8), dimension(:,:),pointer :: bDown3d ! downward component of magnetic field
493 :
494 : real(r8), dimension(pcols,pver) :: sourceR ! R term of source g4 calculation
495 : real(r8), dimension(pcols,pver) :: sourceEff ! Efficiency term of source g4 calculation
496 :
497 21888 : real(r8), dimension(:,:),pointer :: rairvi ! Constituent dependent gas constant
498 :
499 21888 : real(r8), dimension(:,:),pointer :: dipMag ! dip angle for each column (radians)
500 :
501 : real(r8), parameter :: rMassN2 = 28._r8 ! N2 molecular weight kg/kmol
502 : real(r8) :: rMass ! Constituent molecular weight kg/kmol
503 :
504 21888 : real(r8), dimension(:,:),pointer :: mmrN2 ! N2 mass mixing ratio kg/kg
505 : real(r8), dimension(pcols,pver) :: mmrO2 ! O2 mass mixing ratio kg/kg
506 : real(r8), dimension(pcols,pver) :: mmrO1 ! O mass mixing ratio kg/kg
507 :
508 : real(r8), dimension(pcols,pver) :: mmr ! Constituent mass mixing ratio kg/kg
509 :
510 21888 : real(r8), dimension(:,:),pointer :: ndensN2 ! N2 number density (cm-3)
511 : real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3)
512 : real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3)
513 : real(r8), dimension(pcols,pver) :: ndensN1 ! N number density (cm-3)
514 : real(r8), dimension(pcols,pver) :: ndensE ! E electron number density (cm-3)
515 :
516 21888 : real(r8), dimension(:,:) ,pointer :: ndens ! Constituent number density (cm-3)
517 :
518 21888 : real(r8), dimension(:,:,:),pointer :: ionRates ! Pointer to ionization rates for O+,O2+,N+,N2+,NO+ in pbuf (s-1)
519 : ! (from modules mo_jeuv and mo_jshort)
520 :
521 21888 : real(r8), dimension(:,:,:),pointer :: ionPRates ! ionization rates temporary array (s-1 cm-3)
522 21888 : real(r8), dimension(:,:) ,pointer :: sumIonPRates ! Sum of ionization rates for O+,O2+,N+,N2+,NO+ (s-1 cm-3)
523 :
524 21888 : real(r8), dimension(:,:) ,pointer :: aurIPRateSum ! Auroral ion production sum for O2+,O+,N2+
525 : ! (s-1 cm-3 from module mo_aurora)
526 :
527 21888 : real(r8), dimension(:,:) ,pointer :: sourceg4 ! g4 source term for electron/ion temperature update
528 :
529 : real(r8) :: calDay ! current calendar day
530 :
531 : !real(r8), dimension(pcols,pver) :: tempout ! temporary scratch output array
532 :
533 : real(r8), dimension(pcols,pver) :: zGeom ! Geometric altitude (cm)
534 : real(r8), dimension(pcols,pver) :: zThickness ! Geometric altitude thickness (cm)
535 : real(r8), dimension(pcols) :: eColDens ! Electron column density (TECU = 1E16 m-2)
536 :
537 : !--------------------------------------------------------------------------------
538 :
539 21888 : sourceR = 0._r8
540 21888 : sourceEff = 0._r8
541 :
542 21888 : mmrN2 => istate%n2_mmr
543 48394368 : mmrN2 = 0._r8
544 21888 : mmrO2 = 0._r8
545 21888 : mmrO1 = 0._r8
546 21888 : mmr = 0._r8
547 :
548 21888 : ndensO2(:,:) = 0._r8
549 21888 : ndensO1(:,:) = 0._r8
550 21888 : ndensN1(:,:) = 0._r8
551 21888 : ndensE(:,:) = 0._r8
552 :
553 21888 : sourceR(:,:) = 0._r8
554 21888 : sourceEff(:,:) = 0._r8
555 :
556 : !tempout(:,:) = 0._r8
557 :
558 : !--------------------------------------------------------------------------------------
559 : ! Get lchnk from state
560 : !--------------------------------------------------------------------------------------
561 21888 : lchnk = state%lchnk
562 21888 : ncol = state%ncol
563 :
564 : !------------------------------------------------------------------------------------------------------
565 : ! Set the bottom of the ionosphere calculations at around 50 Pascals or 0.5 hectopascals(millibars).
566 : ! teTiBotPres is in Pascals.
567 : !------------------------------------------------------------------------------------------------------
568 21888 : teTiBot = press_lim_idx(teTiBotPres, top=.false.)
569 :
570 : !----------------------------------------------------------------
571 : ! Get latitude and longitude of each column in this chunk
572 : !----------------------------------------------------------------
573 21888 : geoLatR => state%lat(1:ncol)
574 21888 : geoLonR => state%lon(1:ncol)
575 :
576 : !-------------------------------------------------------------------------------------------------------
577 : ! Need to get midpoint and interface pressure and neutral temperature from state structure (pcols,pver)
578 : !-------------------------------------------------------------------------------------------------------
579 21888 : pMid => state%pmid(1:ncol,1:pver)
580 21888 : tN => state%t(1:ncol,1:pver)
581 :
582 21888 : tNInt => istate%tNInt(1:ncol,1:pverp)
583 21888 : cosZenAngR => istate%cosZenAngR(1:ncol)
584 21888 : zenAngD => istate%zenAngD(1:ncol)
585 :
586 21888 : bNorth3d => istate%bNorth3d(1:ncol,1:pver)
587 21888 : bEast3d => istate%bEast3d(1:ncol,1:pver)
588 21888 : bDown3d => istate%bDown3d(1:ncol,1:pver)
589 :
590 21888 : rairvi => istate%rairvi(1:ncol,1:pverp)
591 :
592 21888 : dipMag => istate%dipMag(1:ncol,1:pver)
593 :
594 21888 : ndensN2 => istate%ndensN2(1:ncol,1:pver)
595 :
596 21888 : ionPRates => istate%ionPRates(1:ncol,1:pver,1:nIonRates)
597 21888 : sumIonPRates => istate%sumIonPRates(1:ncol,1:pver)
598 :
599 21888 : sourceg4 => istate%sourceg4(1:ncol,1:pver)
600 :
601 : !-------------------------------------------------------------------------------------
602 : ! Calculate neutral temperature on interface levels. tN vertical dimension is pver
603 : !-------------------------------------------------------------------------------------
604 2845440 : do iVer = 2, pver
605 :
606 36728064 : do iCol = 1, ncol
607 :
608 36706176 : tNInt(iCol,iVer) = 0.5_r8 * tN(iCol,iVer) + 0.5_r8 * tN(iCol,iVer-1)
609 :
610 : enddo
611 : enddo
612 :
613 284544 : do iCol = 1, ncol
614 284544 : tNInt(iCol,1) = 1.5_r8 * tNInt(iCol,2) - 0.5_r8 * tNInt(iCol,3)
615 : enddo
616 284544 : do iCol = 1, ncol
617 284544 : tNInt(iCol,pverp) = 1.5_r8 * tNInt(iCol,pver) - 0.5_r8 * tNInt(iCol,pver-1)
618 : enddo
619 :
620 : !--------------------------------------------------------------
621 : ! Get zenith angle
622 : !--------------------------------------------------------------
623 21888 : calDay = get_curr_calday()
624 21888 : call zenith(calDay,geoLatR(1:ncol),geoLonR(1:ncol),cosZenAngR(1:ncol),ncol)
625 :
626 284544 : do iCol = 1, ncol
627 :
628 284544 : zenAngD(iCol) = ACOS(cosZenAngR(iCol)) * rads2Degs
629 :
630 : enddo
631 :
632 : !---------------------------------------------------------------------------------------
633 : ! Expand magnetic field components in vertical to make 3D, pcols,pver,begchunk:endchunk
634 : ! These are used in calculation of magnetic dip angle and magnetic declination angle so
635 : ! store in local ionosphere module structure.
636 : !---------------------------------------------------------------------------------------
637 2867328 : do iVer = 1, pver
638 :
639 37012608 : do iCol = 1, ncol
640 :
641 34145280 : bNorth3d(iCol,iVer) = bnorth(iCol,lchnk)
642 34145280 : bEast3d(iCol,iVer) = beast(iCol,lchnk)
643 36990720 : bDown3d(iCol,iVer) = bdown(iCol,lchnk)
644 :
645 : enddo
646 :
647 : enddo
648 :
649 : !------------------------------------------------------------------------
650 : ! Get constituent dependent gas constant and derive on interface levels
651 : !------------------------------------------------------------------------
652 2845440 : do iVer = 2, pver
653 36728064 : do iCol = 1, ncol
654 36706176 : rairvi(iCol,iVer) = 0.5_r8 * rairv(iCol,iVer-1,lchnk) + 0.5_r8 * rairv(iCol,iVer,lchnk)
655 : enddo
656 : enddo
657 :
658 284544 : do iCol = 1, ncol
659 284544 : rairvi(iCol,1) = 1.5_r8 * rairvi(iCol,2) - 0.5_r8 * rairvi(iCol,3)
660 : enddo
661 284544 : do iCol = 1, ncol
662 284544 : rairvi(iCol,pverp) = 1.5_r8 * rairvi(iCol,pver) - 0.5_r8 * rairvi(iCol,pver-1)
663 : enddo
664 :
665 : !-------------------------------------------------------------------------------
666 : ! Need to get dip angle from magnetic field components
667 : !-------------------------------------------------------------------------------
668 2867328 : do iVer = 1, pver
669 37012608 : do iCol = 1, ncol
670 34145280 : dipMag(iCol,iVer) = ATAN(bDown3d(iCol,iVer) / SQRT(bNorth3d(iCol,iVer)**2 + bEast3d(iCol,iVer)**2))
671 34145280 : if (dipMag(iCol,iVer) < 0.17_r8 .and. dipMag(iCol,iVer) > 0._r8 ) dipMag(iCol,iVer) = 0.17_r8
672 36990720 : if (dipMag(iCol,iVer) > -0.17_r8 .and. dipMag(iCol,iVer) < 0._r8 ) dipMag(iCol,iVer) = 0.17_r8
673 : enddo
674 : enddo
675 :
676 : !-------------------------------------------------------------------------------------------
677 : ! Set up constituents to be accessed here from pbuf or state%q.
678 : !-------------------------------------------------------------------------------------------
679 218880 : cCnst = (/'O ','O2 ','NO ','H ','N ','e ','Op ','O2p','NOp'/)
680 :
681 218880 : do iCnst = 1, nCnst
682 :
683 : !--------------------------------------
684 : ! Assign density to istate array
685 : !--------------------------------------
686 196992 : if (cCnst(iCnst) == 'O ') ndens => istate%ndensO1(1:ncol,1:pver)
687 196992 : if (cCnst(iCnst) == 'O2 ') ndens => istate%ndensO2(1:ncol,1:pver)
688 196992 : if (cCnst(iCnst) == 'NO ') ndens => istate%ndensNO(1:ncol,1:pver)
689 196992 : if (cCnst(iCnst) == 'N ') ndens => istate%ndensN1(1:ncol,1:pver)
690 196992 : if (cCnst(iCnst) == 'e ') ndens => istate%ndensE(1:ncol,1:pver)
691 196992 : if (cCnst(iCnst) == 'Op ') ndens => istate%ndensOp(1:ncol,1:pver)
692 196992 : if (cCnst(iCnst) == 'O2p') ndens => istate%ndensO2p(1:ncol,1:pver)
693 196992 : if (cCnst(iCnst) == 'NOp') ndens => istate%ndensNOp(1:ncol,1:pver)
694 :
695 : !-------------------------------------------------------------------------------------------
696 : ! Set flag and get field mmr whether each constituent is short-lived(pbuf) or not(state%q).
697 : !-------------------------------------------------------------------------------------------
698 196992 : call cnst_get_ind( TRIM(cCnst(iCnst)), indxCnst, abort=.false. )
699 196992 : if (indxCnst < 0) then
700 87552 : indxSlvd = slvd_index( TRIM(cCnst(iCnst)) )
701 87552 : if (indxSLvd > 0) then
702 350208 : call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxSLvd/), kount=(/pcols,pver,1/) )
703 148050432 : mmr(1:ncol,1:pver) = mmrP(1:ncol,1:pver)
704 87552 : sIndx = get_spc_ndx( TRIM(cCnst(iCnst)) )
705 87552 : rMass = adv_mass(sIndx)
706 : endif
707 : else
708 185063040 : mmr(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxCnst)
709 109440 : rMass = cnst_mw(indxCnst)
710 : endif
711 :
712 : !--------------------------------------------------------------------------------------------------------------
713 : ! Need to get number density (cgs units) from mass mixing ratio. mbarv is kg/mole, same as rMass units
714 : ! kg/kg * (kg/mole)/(kg/mole) * (Pa or N/m*m)/((Joules/K or N*m/K) * (K)) = m-3 * 1E-06 = cm-3
715 : !---------------------------------------------------------------------------------------------------------------
716 787968 : ndens(1:ncol,1:pver) = mmr(1:ncol,1:pver) * mbarv(1:ncol,1:pver,lchnk) / rMass * &
717 666817920 : pMid(1:ncol,1:pver) / (kboltz * tN(1:ncol,1:pver)) * 1.E-06_r8
718 :
719 196992 : if (cCnst(iCnst) == 'O ') then
720 37012608 : mmrO1(1:ncol,1:pver) = mmr(1:ncol,1:pver)
721 37012608 : ndensO1(1:ncol,1:pver) = ndens(1:ncol,1:pver)
722 : endif
723 196992 : if (cCnst(iCnst) == 'O2 ') then
724 37012608 : mmrO2(1:ncol,1:pver) = mmr(1:ncol,1:pver)
725 37012608 : ndensO2(1:ncol,1:pver) = ndens(1:ncol,1:pver)
726 : endif
727 37187712 : if (cCnst(iCnst) == 'N ') ndensN1(1:ncol,1:pver) = ndens(1:ncol,1:pver)
728 37187712 : if (cCnst(iCnst) == 'e ') ndensE(1:ncol,1:pver) = ndens(1:ncol,1:pver)
729 :
730 : !----------------------------------------------------------------------------
731 : ! Calculate N2 density from O2 and O and assign to istate array
732 : !----------------------------------------------------------------------------
733 218880 : if (iCnst == nCnst) then
734 :
735 37012608 : mmrN2(1:ncol,1:pver) = 1._r8 - (mmrO2(1:ncol,1:pver) + mmrO1(1:ncol,1:pver))
736 37012608 : mmrN2(1:ncol,1:pver) = MAX(1.e-20_r8,mmrN2(1:ncol,1:pver))
737 43776 : ndensN2(1:ncol,1:pver) = mmrN2(1:ncol,1:pver) * mbarv(1:ncol,1:pver,lchnk) / rMassN2 * &
738 74047104 : pMid(1:ncol,1:pver) / (kboltz * tN(1:ncol,1:pver)) * 1.E-06_r8
739 :
740 : endif
741 :
742 : enddo ! nCnst
743 :
744 21888 : if (hist_fld_active('ElecColDens')) then
745 : !---------------------------------------
746 : ! Calculate electron column density
747 : !---------------------------------------
748 : !------------------------------------------------------------------------------
749 : ! Convert geopotential altitude in meters to geometric altitude in centimeters
750 : !------------------------------------------------------------------------------
751 37012608 : zGeom(1:ncol,1:pver) = state%zm(1:ncol,1:pver) * (1._r8 + state%zm(1:ncol,1:pver) / rearth) * 100._r8
752 :
753 : !------------------------------------------------------------
754 : ! Calculate vertical thickness at each level in centimeters
755 : !------------------------------------------------------------
756 2823552 : do iVer = 2, pver-1
757 :
758 36443520 : zThickness(1:ncol,iVer) = (zGeom(1:ncol,iVer-1) - zGeom(1:ncol,iVer+1)) / 2._r8
759 :
760 : enddo
761 :
762 284544 : zThickness(1:ncol,1) = (1.5_r8 * zThickness(1:ncol,2)) - (0.5_r8 * zThickness(1:ncol,3))
763 284544 : zThickness(1:ncol,pver) = (1.5_r8 * zThickness(1:ncol,pver-1)) - (0.5_r8 * zThickness(1:ncol,pver-2))
764 :
765 : !----------------------------------------------------------------------------------
766 : ! Calculate electron column density converting from cm-2 to TEC units (1E16 m-2)
767 : ! and make available for history output
768 : !----------------------------------------------------------------------------------
769 34429824 : eColDens(1:ncol) = sum(ndensE(1:ncol,:) * zThickness(1:ncol,:), dim=2) / 1.E12_r8
770 :
771 21888 : call outfld('ElecColDens', eColDens, pcols, lchnk)
772 : endif
773 :
774 : !------------------------------------------------------------------------------------
775 : ! Get ionization rates from physics buffer which were calculated in mo_jeuv and
776 : ! mo_jshort modules. Rates array dimensions are pcols, pver, nIonRates. Units s-1
777 : !------------------------------------------------------------------------------------
778 21888 : indxIR = pbuf_get_index( 'IonRates' )
779 21888 : call pbuf_get_field(pbuf, indxIR, ionRates)
780 :
781 : !----------------------------------------------------------------------------------------------
782 : ! Need to convert these ionization rates to ion production rates by multiplying number density
783 : ! of neutral species appropriate from reactions in mo_jeuv(jeuv) and mo_jshort(jshort)(for NO)
784 : !----------------------------------------------------------------------------------------------
785 2867328 : do iVer = 1, pver
786 37012608 : do iCol = 1, ncol
787 :
788 409743360 : do iIonR = 1, nIonRates
789 409743360 : IF (iIonR <= 3) then
790 102435840 : ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensO1(iCol,iVer)
791 273162240 : else IF (iIonR == 4) then
792 34145280 : ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensN1(iCol,iVer)
793 239016960 : else IF ((iIonR == 5) .OR. (iIonR >= 7 .AND. iIonR <= 9)) then
794 :
795 136581120 : ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensO2(iCol,iVer)
796 102435840 : else IF (iIonR == 6 .OR. iIonR == 10 .OR. iIonR == 11) then
797 102435840 : ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensN2(iCol,iVer)
798 : endif
799 : enddo
800 :
801 : !----------------------------------------------
802 : ! Sum ion production rates all reactions
803 : !----------------------------------------------
804 412588800 : sumIonPRates(iCol,iVer) = SUM(ionPRates(iCol,iVer,1:11))
805 :
806 : enddo
807 : enddo
808 21888 : if (.not.steady_state_ion_elec_temp) then
809 : !-------------------------------------------------------------------------------------------
810 : ! Get aurora ion production rate sum from physics buffer which were calculated in mo_aurora
811 : ! module. Rate array dimensions are pcols, pver. Units s-1 cm-3
812 : !-------------------------------------------------------------------------------------------
813 21888 : indxAIPRS = pbuf_get_index( 'AurIPRateSum' )
814 21888 : call pbuf_get_field(pbuf, indxAIPRS, aurIPRateSum)
815 :
816 : !-------------------------------------------------------------------------------------------------
817 : ! Calculate electron heating rate which is a source in electron/ion temperature derivation
818 : !-------------------------------------------------------------------------------------------------
819 1860480 : do iVer = 1, teTiBot
820 23923584 : do iCol = 1, ncol
821 44126208 : sourceR(iCol,iVer) = LOG( ndensE(iCol,iVer) / (ndensO2(iCol,iVer) + ndensN2(iCol,iVer) + &
822 44126208 : 0.1_r8 * ndensO1(iCol,iVer)) )
823 : sourceEff(iCol,iVer) = EXP( -(12.75_r8 + 6.941_r8 * sourceR(iCol,iVer) + 1.166_r8 * sourceR(iCol,iVer)**2 + &
824 22063104 : 0.08043_r8 * sourceR(iCol,iVer)**3 + 0.001996_r8 * sourceR(iCol,iVer)**4) )
825 :
826 : !-------------------------------------------------------------------------------
827 : ! Calculate g4 source term for electron temperature update
828 : !-------------------------------------------------------------------------------
829 23901696 : sourceg4(iCol,iVer) = (sumIonPRates(iCol,iVer) + aurIPRateSum(iCol,iVer)) * sourceEff(iCol,iVer)
830 :
831 : enddo
832 :
833 : enddo
834 :
835 21888 : call outfld ('SOURCER' , sourceR , pcols, lchnk)
836 21888 : call outfld ('SOURCEEff' , sourceEff , pcols, lchnk)
837 21888 : call outfld ('AURIPRATESUM', aurIPRateSum, pcols, lchnk)
838 :
839 : !----------------------------------------------------------------------------------------------
840 : ! Get Pedersen and Hall Conductivities from physics buffer which were calculated in iondrag
841 : ! module. Conductivity array dimensions are pcols, pver
842 : !-------------------------------------------------------------------------------
843 21888 : indxSP = pbuf_get_index( 'PedConduct' )
844 21888 : indxSH = pbuf_get_index( 'HallConduct' )
845 21888 : call pbuf_get_field(pbuf, indxSP, sigma_ped)
846 21888 : call pbuf_get_field(pbuf, indxSH, sigma_hall)
847 :
848 : endif
849 :
850 21888 : return
851 :
852 43776 : end subroutine update_istate
853 : !
854 : !===============================================================================
855 :
856 21888 : subroutine update_teti(state, dSETendIn, dSETendOut, ztodt, istate, tE, tI, teTiBot)
857 :
858 : !-----------------------------------------------------------------------
859 : ! Routine to compute the electron and ion temperature
860 : !-----------------------------------------------------------------------
861 :
862 21888 : use physconst, only: gravit ! Gravity (m/s2)
863 : use air_composition, only: rairv, mbarv ! Constituent dependent rair and mbar
864 : use mo_apex, only: alatm
865 :
866 : !------------------------------Arguments--------------------------------
867 :
868 : type(physics_state), intent(in), target :: state ! physics state structure
869 : type(ionos_state), intent(in), target :: istate ! ionosphere state structure
870 :
871 : real(r8), dimension(pcols,pver), intent(in) :: dSETendIn ! dry static energy tendency
872 : real(r8), dimension(pcols,pver), intent(out) :: dSETendOut ! dry static energy tendency
873 :
874 : real(r8), intent(in) :: ztodt ! physics time step
875 :
876 : real(r8), dimension(:,:), pointer, intent(inout) :: tE ! Pointer to electron temperature in pbuf (K)
877 : real(r8), dimension(:,:), pointer, intent(inout) :: tI ! Pointer to ion temperature in pbuf (K)
878 :
879 : integer, intent(in) :: teTiBot ! bottom of ionosphere calculations
880 :
881 : !---------------------------Local storage-------------------------------
882 : integer, parameter :: maxIter = 6 ! maximum number of iterations to solve for electron/ion temperature
883 :
884 : integer :: lchnk ! Chunk number
885 : integer :: ncol ! Number of atmospheric columns
886 : integer :: teTiBotP ! bottom of ionosphere calculations plus one more level
887 :
888 : integer :: iVer ! Counter for vertical loops
889 : integer :: iCol ! Counter for column loops
890 : integer :: iter ! Counter for iteration loop
891 :
892 : real(r8), parameter :: Kec1 = 7.5E5_r8 ! c1 constant for calculation of electron conductivity(Ke)
893 : real(r8), parameter :: Kec2 = 3.22E4_r8 ! c2 constant for calculation of electron conductivity(Ke)
894 : real(r8), parameter :: stepweight = 1.0_r8 ! weight of previous and current times step for diagonals
895 : real(r8), parameter :: sToQConv = 6.24E15_r8 ! Conversion from J/kg/s to ev/g/s
896 :
897 : real(r8), parameter :: lossc5 = 1.21E-4_r8 ! c5 constant needed for loss term g3 for electron temperature update
898 : real(r8), parameter :: lossc7 = 3.6E-2_r8 ! c7 constant needed for loss term g3 for electron temperature update
899 : real(r8), parameter :: lossc9 = 5.7E-4_r8 ! c9 constant needed for loss term g3 for electron temperature update
900 : real(r8), parameter :: lossc13 = 7.E-5_r8 ! c13 constant needed for loss term g3 for electron temperature update
901 :
902 : real(r8), parameter :: lossc4pCoef = 1.77E-19_r8
903 : real(r8), parameter :: lossc6pCoef = 1.21E-18_r8
904 : real(r8), parameter :: lossc8pCoef = 7.9E-19_r8
905 : real(r8), parameter :: lossc10pCoef = 1.3E-4_r8
906 : real(r8), parameter :: lossc11pCoef = 3.125E-21_r8
907 : real(r8), parameter :: lossc12pCoef = 3.4E-12_r8
908 : real(r8), parameter :: lossc14pCoef = 1.57E-12_r8
909 : real(r8), parameter :: lossc15pCoef = 2.9E-14_r8
910 : real(r8), parameter :: lossc16pCoef = 6.9E-14_r8
911 : real(r8), parameter :: lossc3pC1 = 3.2E-8_r8
912 : real(r8), parameter :: lossc3pC2 = 15._r8
913 : real(r8), parameter :: lossc3pC3 = 0.53_r8
914 :
915 : real(r8), parameter :: losscinCoef1 = 6.6e-14_r8
916 : real(r8), parameter :: losscinCoef2 = 5.8e-14_r8
917 : real(r8), parameter :: losscinCoef3 = 0.21e-14_r8
918 : real(r8), parameter :: losscinCoef4 = 5.9e-14_r8
919 : real(r8), parameter :: losscinCoef5 = 5.45e-14_r8
920 : real(r8), parameter :: losscinCoef6 = 4.5e-14_r8
921 : real(r8), parameter :: losscinCoef7 = 5.8e-14_r8
922 : real(r8), parameter :: losscinCoef8 = 0.14e-14_r8
923 : real(r8), parameter :: losscinCoef9 = 4.4e-14_r8
924 :
925 : real(r8), parameter :: FeDCoef1 = -9.0E+7_r8
926 : real(r8), parameter :: FeDCoef2 = 4.0E+7_r8
927 :
928 : real(r8), parameter :: losscACoef1 = 5.71E-8_r8
929 : real(r8), parameter :: losscACoef2 = -3352.6_r8
930 : real(r8), parameter :: losscACoef3 = 2.0E-7_r8
931 : real(r8), parameter :: losscACoef4 = -4605.2_r8
932 : real(r8), parameter :: losscACoef5 = 2.53E-6_r8
933 : real(r8), parameter :: losscACoef6 = -17620._r8
934 :
935 : real(r8), parameter :: loss10pCoef = 3200._r8
936 : real(r8), parameter :: lossc12pC1 = 0.4_r8
937 : real(r8), parameter :: lossc12pC2 = 150._r8
938 :
939 : real(r8), parameter :: losscf2dC1 = 2.4E+4_r8
940 : real(r8), parameter :: losscf2dC2 = 0.3_r8
941 : real(r8), parameter :: losscf2dC3 = 1500._r8
942 : real(r8), parameter :: losscf2dC4 = 1.947E-5_r8
943 : real(r8), parameter :: losscf2dC5 = 4000._r8
944 :
945 : real(r8), parameter :: losscf2C1 = 3000._r8
946 :
947 : real(r8), parameter :: losscf3c1 = -22713._r8
948 :
949 : real(r8), parameter :: f1Ted1C1 = 2.82E-17_r8
950 : real(r8), parameter :: f1Ted1C2 = 3.41E-21_r8
951 :
952 : real(r8), parameter :: f1Ted2C1 = 2.2E-16_r8
953 : real(r8), parameter :: f1Ted2C2 = 7.92E-18_r8
954 :
955 : real(r8), parameter :: f1Ted3C1 = 1.1E-16_r8
956 : real(r8), parameter :: f1Ted3C2 = 5.7E-4_r8
957 :
958 : real(r8) :: wrk1 ! 2/3/kboltz_ev
959 : real(r8) :: FeDB ! B term of electron heat flux of UB
960 : real(r8) :: FeD ! Day time flux
961 : real(r8) :: FeN ! Night time flux
962 : real(r8) :: f1Ted1 ! d1 of f1(Te) calculation used to get electron conductivity
963 : real(r8) :: f1Ted2 ! d2 of f1(Te) calculation used to get electron conductivity
964 : real(r8) :: f1Ted3 ! d3 of f1(Te) calculation used to get electron conductivity
965 : real(r8) :: f1Te
966 :
967 21888 : real(r8), dimension(:,:), pointer :: pMid ! Midpoint pressure (Pa)
968 21888 : real(r8), dimension(:,:), pointer :: tN ! Neutral temperature (K)
969 : real(r8), dimension(pcols,pver) :: tEPrevI ! Electron temperature from previous iteration (K)
970 :
971 21888 : real(r8), dimension(:,:), pointer :: pInt ! Interface pressure (Pa)
972 21888 : real(r8), dimension(:,:), pointer :: tNInt ! Interface Temperture (K)
973 21888 : real(r8), dimension(:,:), pointer :: rairvi ! Constituent dependent gas constant on interface levels
974 :
975 21888 : real(r8), dimension(:,:), pointer :: ndensN2 ! N2 number density (cm-3)
976 21888 : real(r8), dimension(:,:), pointer :: ndensO2 ! O2 number density (cm-3)
977 21888 : real(r8), dimension(:,:), pointer :: ndensO1 ! O number density (cm-3)
978 21888 : real(r8), dimension(:,:), pointer :: ndensE ! E electron number density (cm-3)
979 21888 : real(r8), dimension(:,:), pointer :: ndensOp ! O plus number density (cm-3)
980 21888 : real(r8), dimension(:,:), pointer :: ndensO2p ! O2 plus ion number density (cm-3)
981 21888 : real(r8), dimension(:,:), pointer :: ndensNOp ! NO plus ion number density (cm-3)
982 :
983 21888 : real(r8), dimension(:,:), pointer :: sourceg4 ! g4 source term for electron/ion temperature update
984 :
985 21888 : real(r8), dimension(:,:), pointer :: dipMag ! dip angle for each column (radians)
986 : real(r8), dimension(pcols) :: dlatm ! magnetic latitude of each phys column (degrees)
987 :
988 21888 : real(r8), dimension(:), pointer :: zenAngD ! zenith angle (degrees)
989 :
990 : real(r8), dimension(pcols) :: FeUB ! electron heat flux at upper boundary
991 :
992 : real(r8), dimension(pver) :: sqrtTE ! Square root of electron temperature
993 :
994 : real(r8), dimension(pver) :: Ke ! electron conductivity
995 :
996 : real(r8), dimension(pverp) :: Kei ! electron conductivity interface levels
997 :
998 : real(r8), dimension(pcols,pver) :: lossc4p ! c4 prime of Lc(eN2) component of loss term
999 : real(r8), dimension(pcols,pver) :: lossceN2 ! Lc(eN2) component of loss term equation
1000 :
1001 : real(r8), dimension(pcols,pver) :: lossc6p ! c6 prime of Lc(eO2) component of loss term equation
1002 : real(r8), dimension(pcols,pver) :: lossceO2 ! Lc(eO2) component of loss term equation
1003 :
1004 : real(r8), dimension(pcols,pver) :: lossc8p ! c8 prime of Lc(eO) component of loss term equation
1005 : real(r8), dimension(pcols,pver) :: lossceO1 ! Lc(eO) component of loss term equation
1006 :
1007 : real(r8), dimension(pcols,pver) :: lossc10p ! c10 prime of Lc(eN2) component of loss term equation
1008 : real(r8), dimension(pcols,pver) :: losscA ! A of Lc(eN2)v component of loss term equation
1009 : real(r8), dimension(pcols,pver) :: tENDiff ! Difference between electron and neutral temperatures
1010 : real(r8), dimension(pcols,pver) :: lossceN2v ! Lc(eN2)v component of loss term equation
1011 :
1012 : real(r8), dimension(pcols,pver) :: lossc11p ! c11 prime of Lc(eO2)v component of loss term equation
1013 : real(r8), dimension(pcols,pver) :: lossceO2v ! Lc(eO2)v component of loss term equation
1014 :
1015 : real(r8), dimension(pcols,pver) :: lossc12p ! c12 prime of Lc(eO)f component of loss term equation
1016 : real(r8), dimension(pcols,pver) :: lossceOf ! Lc(eO)f component of loss term equation
1017 :
1018 : real(r8), dimension(pcols,pver) :: lossc14p ! c14 prime of Lc(eO)1D component of loss term equation
1019 : real(r8), dimension(pcols,pver) :: losscf2d ! d of f2 of Lc(eO)1D component of loss term equation
1020 : real(r8), dimension(pcols,pver) :: losscf2 ! f2 of Lc(eO)1D component of loss term equation
1021 : real(r8), dimension(pcols,pver) :: losscf3 ! f3 of Lc(eO)1D component of loss term equation
1022 : real(r8), dimension(pcols,pver) :: lossceO1D ! Lc(eO)1D component of loss term equation
1023 :
1024 : real(r8), dimension(pcols,pver) :: lossc15p ! c15 prime of Lc(eN2)Rot component of loss term equation
1025 : real(r8), dimension(pcols,pver) :: lossceN2Rot ! Lc(eN2)Rot component of loss term equation
1026 :
1027 : real(r8), dimension(pcols,pver) :: lossc16p ! c16 prime of Lc(eO2)Rot component of loss term equation
1028 : real(r8), dimension(pcols,pver) :: lossceO2Rot ! Lc(eO2)Rot component of loss term equation
1029 :
1030 : real(r8), dimension(pcols,pver) :: lossc3p ! c3 prime of Lc(ei) component of loss term equation
1031 : real(r8), dimension(pcols,pver) :: losscei ! Lc(ei) component of loss term equation
1032 : real(r8), dimension(pcols,pver) :: losscin ! ion-neutral heating coeff.
1033 :
1034 : real(r8), dimension(pcols,pver) :: lossg3 ! g3 loss term for Te tendency
1035 :
1036 : real(r8), dimension(pcols,pverp) :: delZi ! Delta z: interfaces
1037 : real(r8), dimension(pcols,pver) :: delZ ! Delta z: midpoints
1038 :
1039 : real(r8), dimension(pcols,pver) :: qjoule ! joule heating
1040 : real(r8), dimension(pcols,pver) :: qen ! electron-neutral heating (units: ev/g/s)
1041 : real(r8), dimension(pcols,pver) :: qei ! electron-ion Coulomb heating (units: ev/g/s)
1042 : real(r8), dimension(pcols,pver) :: qin ! ion-neutral heating (units: ev/g/s)
1043 : real(r8), dimension(pcols,pver) :: rho ! mass density
1044 :
1045 : real(r8), dimension(pcols,pver) :: wrk2
1046 :
1047 43776 : real(r8), dimension(teTiBot) :: subdiag ! subdiagonal values for Te tendency solving
1048 43776 : real(r8), dimension(teTiBot) :: superdiag ! superdiagonal values for Te tendency solving
1049 43776 : real(r8), dimension(teTiBot) :: diag ! diagonal values for Te tendency solving
1050 43776 : real(r8), dimension(teTiBot) :: rHS ! RHS of electron temperature update
1051 43776 : real(r8), dimension(teTiBot) :: tETemp ! temporary electron temperature array for input to tridag
1052 :
1053 : logical, dimension(pcols) :: colConv ! flag for column converging
1054 : logical :: converged ! Flag for convergence in electron temperature
1055 : ! calculation iteration loop
1056 : real(r8) :: qrate(pcols,pver) ! heating rate diagnostic
1057 :
1058 : !---------------------------------------------------------------------------------------------------------
1059 : ! Initialize arrays to zero and column convergence logical to .false.
1060 : !---------------------------------------------------------------------------------------------------------
1061 :
1062 21888 : sqrtTE(:) = 0._r8
1063 21888 : Ke(:) = 0._r8
1064 21888 : Kei(:) = 0._r8
1065 21888 : lossc4p(:,:) = 0._r8
1066 21888 : lossceN2(:,:) = 0._r8
1067 21888 : lossc6p(:,:) = 0._r8
1068 21888 : lossceO2(:,:) = 0._r8
1069 21888 : lossc8p(:,:) = 0._r8
1070 21888 : lossceO1(:,:) = 0._r8
1071 21888 : lossc10p(:,:) = 0._r8
1072 21888 : losscA(:,:) = 0._r8
1073 21888 : tENDiff(:,:) = 0._r8
1074 21888 : lossceN2v(:,:) = 0._r8
1075 21888 : lossc11p(:,:) = 0._r8
1076 21888 : lossceO2v(:,:) = 0._r8
1077 21888 : lossc12p(:,:) = 0._r8
1078 21888 : lossceOf(:,:) = 0._r8
1079 21888 : lossc14p(:,:) = 0._r8
1080 21888 : losscf2d(:,:) = 0._r8
1081 21888 : losscf2(:,:) = 0._r8
1082 21888 : losscf3(:,:) = 0._r8
1083 21888 : lossceO1D(:,:) = 0._r8
1084 21888 : lossc15p(:,:) = 0._r8
1085 21888 : lossceN2Rot(:,:) = 0._r8
1086 21888 : lossc16p(:,:) = 0._r8
1087 21888 : lossceO2Rot(:,:) = 0._r8
1088 21888 : lossc3p(:,:) = 0._r8
1089 21888 : losscei(:,:) = 0._r8
1090 21888 : losscin(:,:) = 0._r8
1091 21888 : lossg3(:,:) = 0._r8
1092 21888 : delZi(:,:) = 0._r8
1093 21888 : delZ(:,:) = 0._r8
1094 1860480 : subDiag(:) = 0._r8
1095 1860480 : superDiag(:) = 0._r8
1096 1860480 : diag(:) = 0._r8
1097 1860480 : rHS(:) = 0._r8
1098 1860480 : teTemp(:) = 0._r8
1099 21888 : qjoule(:,:) = 0._r8
1100 21888 : qei(:,:) = 0._r8
1101 21888 : qen(:,:) = 0._r8
1102 21888 : qin(:,:) = 0._r8
1103 21888 : rho(:,:) = 0._r8
1104 21888 : dSETendOut = 0._r8
1105 21888 : colConv(:) = .false.
1106 :
1107 : !--------------------------------------------------------------------------------------
1108 : ! Get lchnk and ncol from state
1109 : !--------------------------------------------------------------------------------------
1110 21888 : lchnk = state%lchnk
1111 21888 : ncol = state%ncol
1112 :
1113 : !-------------------------------------------
1114 : ! Calculate some commonly used variables
1115 : !-------------------------------------------
1116 21888 : wrk1 = 2._r8 / 3._r8/ kboltz_ev
1117 21888 : teTiBotP = teTiBot + 1
1118 :
1119 : !-------------------------------------------------------------------------------------------------------
1120 : ! Need to get midpoint and interface pressure and neutral temperature from state structure (ncol,teTiBot)
1121 : !-------------------------------------------------------------------------------------------------------
1122 21888 : pMid => state%pmid(1:ncol,1:pver)
1123 21888 : tN => state%t(1:ncol,1:pver)
1124 37012608 : rho(1:ncol,1:pver) = pMid(1:ncol,1:pver)/rairv(1:ncol,1:pver,lchnk)/tN(1:ncol,1:pver) * 1.E-3_r8 ! convert to g/cm3
1125 :
1126 23923584 : qjoule(1:ncol,1:teTiBot) = dSETendIn(1:ncol,1:teTiBot) * sToQConv ! convert from J/kg/s to ev/g/s
1127 :
1128 21888 : pInt => state%pint(1:ncol,1:pverp)
1129 21888 : tNInt => istate%tNInt(1:ncol,1:pverp)
1130 21888 : rairvi => istate%rairvi(1:ncol,1:pverp)
1131 :
1132 : !----------------------------------------------------------------
1133 : ! Get variables needed from the ionosphere state structure
1134 : !----------------------------------------------------------------
1135 21888 : ndensO2 => istate%ndensO2(1:ncol,1:pver)
1136 21888 : ndensO1 => istate%ndensO1(1:ncol,1:pver)
1137 21888 : ndensE => istate%ndensE(1:ncol,1:pver)
1138 21888 : ndensOp => istate%ndensOp(1:ncol,1:pver)
1139 21888 : ndensO2p => istate%ndensO2p(1:ncol,1:pver)
1140 21888 : ndensNOp => istate%ndensNOp(1:ncol,1:pver)
1141 21888 : ndensN2 => istate%ndensN2(1:ncol,1:pver)
1142 :
1143 21888 : sourceg4 => istate%sourceg4(1:ncol,1:pver)
1144 :
1145 21888 : dipMag => istate%dipMag(1:ncol,1:pver)
1146 :
1147 21888 : zenAngD => istate%zenAngD(1:ncol)
1148 :
1149 : !-------------------------------------------------------------------------------------------------------------------
1150 : ! Set electron temperature limits
1151 : !-------------------------------------------------------------------------------------------------------------------
1152 74003328 : tE(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),tE(1:ncol,1:pver))
1153 37012608 : tE(1:ncol,1:pver) = MIN(temax,tE(1:ncol,1:pver))
1154 :
1155 74003328 : tI(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),ti(1:ncol,1:pver))
1156 74003328 : tI(1:ncol,1:pver) = MIN(ti(1:ncol,1:pver),tE(1:ncol,1:pver))
1157 :
1158 : ! set Te and Ti to Tn below the levels where this module applies
1159 26199936 : tE(1:ncol,teTiBotP:pver) = tN(1:ncol,teTiBotP:pver)
1160 26199936 : tI(1:ncol,teTiBotP:pver) = tN(1:ncol,teTiBotP:pver)
1161 :
1162 23923584 : wrk2(1:ncol,1:teTiBot) = ndensE(1:ncol,1:teTiBot)/wrk1/(SIN(dipMag(1:ncol,1:teTiBot)))**2._r8
1163 :
1164 284544 : dlatm(:ncol) = rads2Degs* alatm(:ncol,lchnk)
1165 :
1166 : !-----------------------------------------------------------------------------
1167 : ! Get terms needed for loss term g3 for electron temperature update which do
1168 : ! not need to be updated in iteration loop.
1169 : !-----------------------------------------------------------------------------
1170 284544 : do iCol = 1, ncol
1171 :
1172 284544 : if (.not. colConv(iCol)) then
1173 22325760 : do iVer = 1, teTiBot
1174 :
1175 22063104 : lossc4p(iCol,iVer) = lossc4pCoef * ndensN2(iCol,iVer) * ndensE(iCol,iVer) ! e-N2 elastic collision
1176 22063104 : lossc6p(iCol,iVer) = lossc6pCoef * ndensO2(iCol,iVer) * ndensE(iCol,iVer) ! e-O2 elastic collision
1177 22063104 : lossc8p(iCol,iVer) = lossc8pCoef * ndensO1(iCol,iVer) * ndensE(iCol,iVer) ! e-O elastic collision
1178 22063104 : lossc10p(iCol,iVer) = lossc10pCoef * ndensN2(iCol,iVer) * ndensE(iCol,iVer) ! e-N2(vib)
1179 22063104 : lossc11p(iCol,iVer) = lossc11pCoef * ndensO2(iCol,iVer) * ndensE(iCol,iVer) ! e-O2(vib)
1180 22063104 : lossc12p(iCol,iVer) = lossc12pCoef * ndensO1(iCol,iVer) * ndensE(iCol,iVer) ! e-O (fine)
1181 22063104 : lossc14p(iCol,iVer) = lossc14pCoef * ndensO1(iCol,iVer) * ndensE(iCol,iVer) ! e-O(1D)
1182 22063104 : lossc15p(iCol,iVer) = lossc15pCoef * ndensN2(iCol,iVer) * ndensE(iCol,iVer) ! e-N2(rot)
1183 22063104 : lossc16p(iCol,iVer) = lossc16pCoef * ndensO2(iCol,iVer) * ndensE(iCol,iVer) ! e-O2(rot)
1184 : lossc3p(iCol,iVer) = lossc3pC1 * lossc3pC2 * (ndensOP(iCol,iVer) + &
1185 22063104 : 0.5_r8 * ndensO2P(iCol,iVer) + lossc3pC3 * ndensNOP(iCol,iVer)) * ndensE(iCol,iVer) ! e-i
1186 :
1187 : losscin(iCol,iVer) = (losscinCoef1*ndensN2(iCol,iVer) + losscinCoef2*ndensO2(iCol,iVer) &
1188 : + losscinCoef3*ndensO1(iCol,iVer)*SQRT(2._r8*tN(iCol,iVer)))*ndensOP(iCol,iVer) &
1189 : +(losscinCoef4*ndensN2(iCol,iVer) + losscinCoef5*ndensO2(iCol,iVer) &
1190 : + losscinCoef6*ndensO1(iCol,iVer))*ndensNOP(iCol,iVer) &
1191 : +(losscinCoef7*ndensN2(iCol,iVer) + losscinCoef8*ndensO2(iCol,iVer)*SQRT(tN(iCol,iVer)) &
1192 22325760 : + losscinCoef9*ndensO1(iCol,iVer)) * ndensO2P(iCol,iVer)
1193 :
1194 :
1195 : enddo !iVer loop
1196 :
1197 : !----------------------------------------------------------------------------------
1198 : ! Calculate upper boundary heat flux
1199 : !----------------------------------------------------------------------------------
1200 262656 : if (ABS(dlatm(iCol)) < 10.0_r8) then
1201 : FeDB = 0._r8
1202 234555 : else if (ABS(dlatm(iCol)) >= 10.0_r8 .and. ABS(dlatm(iCol)) < 40.0_r8) then
1203 87400 : FeDB = 0.5_r8 * (1._r8 + COS(pi * (ABS(dlatm(iCol)) - 40.0_r8) /30.0_r8))
1204 : else
1205 : FeDB = 1._r8
1206 : end if
1207 :
1208 262656 : FeD = FeDCoef1 * f107 * FeDB
1209 262656 : FeN = .5_r8 * FeD
1210 : !---------------------------------------------------
1211 : ! Set upper boundary condition for right hand side
1212 : !---------------------------------------------------
1213 262656 : if (zenAngD(iCol) <= 80.0_r8) FeUB(iCol) = FeD
1214 262656 : if (zenAngD(iCol) > 80.0_r8 .AND. zenAngD(iCol) < 100.0_r8) FeUB(iCol) = 0.5_r8 * (FeD + FeN) &
1215 : + 0.5_r8 * (FeD - FeN) * &
1216 44172 : COS(pi * ((zenAngD(iCol) - 80.0_r8) / 20.0_r8))
1217 262656 : if (zenAngD(iCol) >= 100.0_r8) FeUB(iCol) = FeN
1218 :
1219 : !------------------------------------------------------------------------------------------
1220 : ! Calculate thickness terms for vertical derivative
1221 : !------------------------------------------------------------------------------------------
1222 22325760 : do iVer = 1, teTiBot
1223 :
1224 44126208 : delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rairv(iCol,iVer,lchnk) * &
1225 66451968 : tN(iCol,iVer) / pMid(iCol,iVer) / gravit
1226 :
1227 : enddo
1228 :
1229 22325760 : do iVer = 2, teTiBotP ! Assuming teTiBotP < pverp
1230 44126208 : delZi(iCol,iVer) = (pMid(iCol,iVer) - pMid(iCol,iVer-1)) * rairvi(iCol,iVer) * &
1231 66451968 : tNInt(iCol,iVer) / pInt(iCol,iVer) / gravit
1232 : enddo
1233 262656 : delZi(iCol,1) = 1.5_r8*delZi(iCol,2) - .5_r8*delZi(iCol,3)
1234 :
1235 : !----------------------------------------------------------
1236 : ! Convert delZ variables from meters to centimeters
1237 : !----------------------------------------------------------
1238 22588416 : delZi(iCol,1:teTiBotP) = delZi(iCol,1:teTiBotP)*100._r8
1239 22325760 : delZ(iCol,1:teTiBot) = delZ(iCol,1:teTiBot)*100._r8
1240 :
1241 : endif ! Column not converged
1242 :
1243 : enddo !iCol loop
1244 :
1245 : !-------------------------------------------------------------------------------------------------------
1246 : ! Iterate to calculate new electron temperature.
1247 : ! Time splitting is used: first solve the heating/cooling equation, then solve the diffusion equations.
1248 : ! Also, set convergence flag to false and iterate until true or 6 iterations, whichever comes first
1249 : !-------------------------------------------------------------------------------------------------------
1250 : converged = .false.
1251 : iter = 0
1252 48746 : do while (.not. converged .and. iter < maxIter)
1253 :
1254 : !--------------------------------------------------------------------------------------------------------
1255 : ! Increment iteration loop counter and save electron temperature from previous iteration for convergence
1256 : ! test at end of interation loop. Also, take square root of electron temperature to be used later
1257 : !--------------------------------------------------------------------------------------------------------
1258 26858 : iter = iter + 1
1259 :
1260 29355794 : tEPrevI(1:ncol,1:teTiBot) = tE(1:ncol,1:teTiBot)
1261 :
1262 : !--------------------------------------------------------------------------------------------------------
1263 : ! Loop over columns then vertical levels and call tridiagonal solver for each column to get electron
1264 : ! temperature
1265 : !--------------------------------------------------------------------------------------------------------
1266 349154 : do iCol = 1, ncol
1267 :
1268 349154 : if (.not. colConv(iCol)) then
1269 :
1270 23076055 : sqrtTE(1:teTiBot) = SQRT(tE(iCol,1:teTiBot))
1271 :
1272 23076055 : do iVer = 1, teTiBot
1273 :
1274 : !-----------------------------------------------------------------------------
1275 : ! Get loss term g3 for electron temperature update. Need to calculate
1276 : ! constituent dependent loss terms which make up g3
1277 : !-----------------------------------------------------------------------------
1278 22804572 : lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer)
1279 22804572 : lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iVer)) * sqrtTE(iVer)
1280 22804572 : lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iVer)
1281 :
1282 22804572 : if (tE(iCol,iVer) < 1000.0_r8) then
1283 13353163 : losscA(iCol,iVer) = losscACoef1 * EXP(losscACoef2 / tE(iCol,iVer))
1284 : endif
1285 22804572 : if (tE(iCol,iVer) >= 1000.0_r8 .AND. tE(iCol,iVer) <= 2000.0_r8) then
1286 5700289 : losscA(iCol,iVer) = losscACoef3 * EXP(losscACoef4 / tE(iCol,iVer))
1287 : endif
1288 22804572 : if (tE(iCol,iVer) > 2000.0_r8) then
1289 3751120 : losscA(iCol,iVer) = losscACoef5 * sqrtTE(iVer) * EXP(losscACoef6 / tE(iCol,iVer))
1290 : endif
1291 :
1292 22804572 : tENDiff(iCol,iVer) = tE(iCol,iVer) - tN(iCol,iVer)
1293 22804572 : if (ABS(tENDiff(iCol,iVer)) < 0.1_r8) tENDiff(iCol,iVer) = 0.1_r8
1294 :
1295 : lossceN2v(iCol,iVer) = lossc10p(iCol,iVer) * losscA(iCol,iVer) * &
1296 22804572 : (1._r8 - EXP(loss10pCoef * (1._r8 / tE(iCol,iVer) - 1._r8 / tN(iCol,iVer)))) / tENDiff(iCol,iVer)
1297 22804572 : lossceO2v(iCol,iVer) = lossc11p(iCol,iVer) * tE(iCol,iVer) * tE(iCol,iVer)
1298 : lossceOf(iCol,iVer) = lossc12p(iCol,iVer) * (1._r8 - lossc13 * tE(iCol,iVer)) * &
1299 22804572 : (lossc12pC1 + lossc12pC2 / tE(iCol,iVer)) / tN(iCol,iVer)
1300 : losscf2d(iCol,iVer) = losscf2dC1 + losscf2dC2 * (tE(iCol,iVer) - losscf2dC3) - &
1301 22804572 : losscf2dC4 * (tE(iCol,iVer) - losscf2dC3) * (tE(iCol,iVer) - losscf2dC5)
1302 22804572 : losscf2(iCol,iVer) = losscf2d(iCol,iVer) * (1._r8 / losscf2C1 - 1._r8 / tE(iCol,iVer))
1303 22804572 : losscf3(iCol,iVer) = losscf3c1 * (1._r8 / tN(iCol,iVer) - 1._r8 / tE(iCol,iVer))
1304 : lossceO1D(iCol,iVer) = lossc14p(iCol,iVer) * EXP(losscf2(iCol,iVer)) * &
1305 22804572 : (1._r8 - EXP(losscf3(iCol,iVer))) / tENDiff(iCol,iVer)
1306 22804572 : lossceN2Rot(iCol,iVer) = lossc15p(iCol,iVer) / sqrtTE(iVer)
1307 22804572 : lossceO2Rot(iCol,iVer) = lossc16p(iCol,iVer) / sqrtTE(iVer)
1308 :
1309 22804572 : losscei(iCol,iVer) = lossc3p(iCol,iVer) / tE(iCol,iVer)**1.5_r8
1310 :
1311 : !------------------------------------------------
1312 : ! Loss term: lossg3*tE/sin(I)^2
1313 : !------------------------------------------------
1314 : lossg3(iCol,iVer) = lossceN2(iCol,iVer) + lossceO2(iCol,iVer) + lossceO1(iCol,iVer) + lossceN2v(iCol,iVer) &
1315 : + lossceO2v(iCol,iVer) + lossceOf(iCol,iVer) + lossceO1D(iCol,iVer) &
1316 23076055 : + lossceN2Rot(iCol,iVer) + lossceO2Rot(iCol,iVer)
1317 :
1318 : enddo !iVer loop
1319 :
1320 : endif ! Column not converged
1321 :
1322 : enddo ! End of column loop
1323 :
1324 : !-----------------------------------------------------
1325 : ! Calculate thermal conductivity of electron gas
1326 : !-----------------------------------------------------
1327 349154 : do iCol = 1, ncol
1328 :
1329 349154 : if (.not. colConv(iCol)) then
1330 :
1331 23076055 : sqrtTE(1:teTiBot) = SQRT(tE(iCol,1:teTiBot))
1332 :
1333 23076055 : do iVer = 1, teTiBot
1334 :
1335 22804572 : f1Ted1 = f1Ted1C1 * sqrtTE(iVer) - f1Ted1C2 * tE(iCol,iVer)**1.5_r8
1336 22804572 : f1Ted2 = f1Ted2C1 + f1Ted2C2 * sqrtTE(iVer)
1337 22804572 : f1Ted3 = f1Ted3C1 * (1._r8 + f1Ted3C2 * tE(iCol,iVer))
1338 :
1339 22804572 : f1Te = ndensN2(iCol,iVer) / ndensE(iCol,iVer) * f1Ted1 + ndensO2(iCol,iVer) / &
1340 22804572 : ndensE(iCol,iVer) * f1Ted2 + ndensO1(iCol,iVer) / ndensE(iCol,iVer) * f1Ted3
1341 :
1342 : !-----------------------------------------------------------------------------
1343 : ! Calculate electron conductivity using parameters set in module and f1(Te)
1344 : !-----------------------------------------------------------------------------
1345 23076055 : Ke(iVer) = Kec1 * tE(iCol,iVer)**2.5_r8 / (1._r8 + Kec2 * tE(iCol,iVer)**2._r8 * f1Te)
1346 :
1347 : enddo !iVer loop
1348 :
1349 : !----------------------------------------------------------------------
1350 : ! Get electron conductivity at interface levels to be used later
1351 : !----------------------------------------------------------------------
1352 22804572 : do iVer = 2,teTiBot
1353 22804572 : Kei(iVer) = SQRT(Ke(iVer-1)*Ke(iVer))
1354 : enddo
1355 271483 : Kei(1) = 1.5_r8*Kei(2)-.5_r8*Kei(3)
1356 271483 : Kei(teTiBotP) = 1.5_r8*Kei(teTiBot)-.5_r8*Kei(teTiBot-1)
1357 :
1358 : !------------------------------------------------------------------------------------------------------
1359 : ! Derive subdiagonal, superdiagonal, and diagonal as input to solver for electron temperature tendency
1360 : !------------------------------------------------------------------------------------------------------
1361 22533089 : do iVer = 2, teTiBot-1
1362 22261606 : subDiag(iVer) = -Kei(iVer) / delZi(iCol,iVer) / delZ(iCol,iVer)
1363 22261606 : superDiag(iVer) = -Kei(iVer+1) / delZi(iCol,iVer+1) / delZ(iCol,iVer)
1364 22261606 : diag(iVer) = wrk2(iCol,iVer)/ztodt + (lossg3(iCol,iVer)+losscei(iCol,iVer))/SIN(dipMag(iCol,iVer))**2._r8 &
1365 22261606 : -subDiag(iVer)-superDiag(iVer)
1366 0 : rHS(iVer) = tE(iCol,iVer) * wrk2(iCol,iVer)/ztodt + sourceg4(iCol,iVer)/SIN(dipMag(iCol,iVer))**2._r8 &
1367 22533089 : +(lossg3(iCol,iVer)*tN(iCol,iVer)+losscei(iCol,iVer)*ti(iCol,iVer))/SIN(dipMag(iCol,iVer))**2._r8
1368 : enddo !iVer loop
1369 :
1370 : !-------------------------------------------------------------------------------------
1371 : ! Calculate diagonal, superdiagonal, and right hand side upper boundary values
1372 : !-------------------------------------------------------------------------------------
1373 271483 : superDiag(1) = -Kei(2) / delZi(iCol,2) / delZ(iCol,1)
1374 271483 : diag(1) = wrk2(iCol,1)/ztodt - superDiag(1)
1375 271483 : rHS(1) = tE(iCol,1) * wrk2(iCol,1) / ztodt - FeUB(iCol) / delZ(iCol,1)
1376 :
1377 : !---------------------------------------------------------------------------------------------
1378 : ! Calculate subdiagonal, diagonal, superdiagonal, and right hand side lower boundary values
1379 : !---------------------------------------------------------------------------------------------
1380 271483 : subDiag(teTiBot) = -Kei(teTiBot) / delZi(iCol,teTiBot) / delZ(iCol,teTiBot)
1381 271483 : superDiag(teTiBot) = -Kei(teTiBotP) / delZi(iCol,teTiBotP) / delZ(iCol,teTiBot)
1382 : diag(teTiBot) = wrk2(iCol,teTiBot)/ztodt &
1383 271483 : + (lossg3(iCol,teTiBot)+losscei(iCol,teTiBot))/SIN(dipMag(iCol,teTiBot))**2._r8 &
1384 271483 : - subDiag(teTiBot)-superDiag(teTiBot)
1385 0 : rHS(teTiBot) = tE(iCol,teTiBot) * wrk2(iCol,teTiBot)/ztodt &
1386 : + sourceg4(iCol,teTiBot)/SIN(dipMag(iCol,teTiBot))**2._r8 &
1387 0 : + (lossg3(iCol,teTiBot)*tN(iCol,teTiBot)+losscei(iCol,teTiBot)*ti(iCol,teTiBot)) &
1388 271483 : / SIN(dipMag(iCol,teTiBot))**2._r8 - superDiag(teTiBot) * tN(iCol,teTiBotP)
1389 :
1390 : !-------------------------------------------------
1391 : ! Call solver to get electron temperature update
1392 : !-------------------------------------------------
1393 271483 : call tridag(subDiag,diag,superDiag,rHS,tETemp,teTiBot)
1394 :
1395 23076055 : tE(iCol,1:teTiBot) = tETemp(1:teTiBot)
1396 23076055 : do iVer = 1,teTiBot
1397 22804572 : tE(iCol,iVer) = min(temax,tE(iCol,iVer))
1398 23076055 : tE(iCol,iVer) = max(tN(iCol,iVer),tE(iCol,iVer))
1399 : enddo
1400 :
1401 : !---------------------------------------------------------------------------------------------------------
1402 : ! Calculate ion temperature from electron temperature, ion-neutral and electron-ion loss terms, neutral
1403 : ! temperature, mass density and joule heating. Set minimum value to neutral temperature and maximum
1404 : ! value to electron temperature for each column and vertical level
1405 : !---------------------------------------------------------------------------------------------------------
1406 23076055 : do iVer = 1,teTiBot
1407 45609144 : ti(iCol,iVer) = (losscei(iCol,iVer) * tE(iCol,iVer) + losscin(iCol,iVer) * tN(iCol,iVer) + &
1408 0 : rho(iCol,iVer) * qjoule(iCol,iVer) * mbarv(iCol,iVer,lchnk) / (mbarv(iCol,iVer,lchnk)+rMassOp)) / &
1409 68413716 : (losscei(iCol,iVer) + losscin(iCol,iVer))
1410 22804572 : ti(iCol,iVer) = max(tN(iCol,iVer),ti(iCol,iVer))
1411 23076055 : ti(iCol,iVer) = min(tE(iCol,iVer),ti(iCol,iVer))
1412 : enddo
1413 :
1414 : !--------------------------------------------------------------------------------------------------------
1415 : ! Check for convergence which is a change of electron temperature ratio to previous loop for all levels
1416 : ! and columns of less than 0.05K. Had to modify this to do convergence check on each column since
1417 : ! checking all columns in a chunk gives different answers depending on number of tasks and tasks per node.
1418 : !--------------------------------------------------------------------------------------------------------
1419 22393028 : if (ALL(ABS(tE(iCol,1:teTiBot) / tEPrevI(iCol,1:teTiBot) - 1._r8) < 0.05_r8)) then
1420 :
1421 262656 : colConv(iCol) = .true.
1422 :
1423 : endif
1424 :
1425 : endif ! Column not converged
1426 :
1427 : enddo ! iCol loop
1428 : !--------------------------------------------------------------
1429 : ! Check to see if all columns have converged and set flag
1430 : !--------------------------------------------------------------
1431 349368 : if (ALL(colConv(1:ncol))) converged = .true.
1432 :
1433 : enddo ! End of iteration loop
1434 :
1435 : !--------------------------------------------------------------------------------------------------------
1436 : ! Calculate electron-neutral heating and electron-ion Coulomb heating. Then update dry static energy.
1437 : !--------------------------------------------------------------------------------------------------------
1438 1860480 : do iVer = 1, teTiBot
1439 23923584 : do iCol = 1, ncol
1440 22063104 : sqrtTE(iVer) = SQRT(tE(iCol,iVer))
1441 22063104 : lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer)
1442 22063104 : lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iVer)) * sqrtTE(iVer)
1443 23901696 : lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iVer)
1444 : enddo
1445 : enddo
1446 :
1447 : qen(1:ncol,1:teTiBot) = (lossceN2(1:ncol,1:teTiBot)+lossceO2(1:ncol,1:teTiBot)+lossceO1(1:ncol,1:teTiBot)+ &
1448 : lossceN2v(1:ncol,1:teTiBot)+lossceO2v(1:ncol,1:teTiBot)+lossceOf(1:ncol,1:teTiBot)+ &
1449 : lossceO1D(1:ncol,1:teTiBot)+lossceN2Rot(1:ncol,1:teTiBot)+lossceO2Rot(1:ncol,1:teTiBot))* &
1450 23923584 : (tE(1:ncol,1:teTiBot)-tN(1:ncol,1:teTiBot)) / rho(1:ncol,1:teTiBot)
1451 23923584 : qei(1:ncol,1:teTiBot) = losscei(1:ncol,1:teTiBot) * (tE(1:ncol,1:teTiBot)-ti(1:ncol,1:teTiBot)) / rho(1:ncol,1:teTiBot)
1452 23923584 : qin(1:ncol,1:teTiBot) = losscin(1:ncol,1:teTiBot) * (tI(1:ncol,1:teTiBot)-tN(1:ncol,1:teTiBot)) / rho(1:ncol,1:teTiBot)
1453 :
1454 23923584 : dSETendOut(1:ncol,1:teTiBot) = (qei(1:ncol,1:teTiBot)+qen(1:ncol,1:teTiBot)) / sToQConv ! J/kg/s
1455 :
1456 37012608 : qrate(:ncol,:) = qen(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
1457 21888 : call outfld ('QEN', qrate, pcols, lchnk)
1458 :
1459 37012608 : qrate(:ncol,:) = qei(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
1460 21888 : call outfld ('QEI', qrate, pcols, lchnk)
1461 :
1462 37012608 : qrate(:ncol,:) = qin(:ncol,:)/sToQConv/cpairv(:ncol,:,lchnk) ! K/s
1463 21888 : call outfld ('QIN', qrate, pcols, lchnk)
1464 :
1465 21888 : return
1466 :
1467 21888 : end subroutine update_teti
1468 :
1469 : !-----------------------------------------------------------------------
1470 : ! Simple tridiagonal solver routine
1471 : !-----------------------------------------------------------------------
1472 :
1473 271483 : SUBROUTINE tridag(a,b,c,r,u,n)
1474 :
1475 : INTEGER,INTENT(IN) :: n
1476 : REAL(r8),INTENT(IN) :: a(n),b(n),c(n),r(n)
1477 : REAL(r8),INTENT(INOUT) :: u(n)
1478 : !------------------------------
1479 : ! Local variables
1480 : !------------------------------
1481 : INTEGER j
1482 542966 : REAL(r8) :: bet,gam(n)
1483 :
1484 271483 : if(b(1).eq.0._r8) call endrun('ion_electron_temp: bt(1)=0 in tridag')
1485 271483 : bet=b(1)
1486 271483 : u(1)=r(1)/bet
1487 22804572 : do j=2,n
1488 22533089 : gam(j)=c(j-1)/bet
1489 22533089 : bet=b(j)-a(j)*gam(j)
1490 22533089 : if(bet.eq.0._r8) call endrun('ion_electron_temp: bet=0 in tridag')
1491 22804572 : u(j)=(r(j)-a(j)*u(j-1))/bet
1492 : end do
1493 :
1494 22804572 : do j=n-1,1,-1
1495 22804572 : u(j)=u(j)-gam(j+1)*u(j+1)
1496 : end do
1497 :
1498 271483 : return
1499 :
1500 : END SUBROUTINE tridag
1501 :
1502 : end module ion_electron_temp
|
