Line data Source code
1 : module ionosphere_interface
2 :
3 : use shr_kind_mod, only: r8 => shr_kind_r8, cl=>shr_kind_cl
4 : use cam_abortutils, only: endrun
5 : use ppgrid, only: begchunk, endchunk, pcols, pver
6 : use phys_grid, only: get_ncols_p
7 :
8 : use dpie_coupling, only: d_pie_init
9 : use dpie_coupling, only: d_pie_epotent
10 : use dpie_coupling, only: d_pie_coupling ! WACCM-X ionosphere/electrodynamics coupling
11 : use short_lived_species, only: slvd_index, slvd_pbf_ndx => pbf_idx ! Routines to access short lived species
12 :
13 : use chem_mods, only: adv_mass ! Array holding mass values for short lived species
14 : use mo_chem_utls, only: get_spc_ndx ! Routine to get index of adv_mass array for short lived species
15 : use physics_buffer, only: pbuf_get_chunk, pbuf_get_field
16 : use physics_buffer, only: pbuf_get_index
17 :
18 : use constituents, only: cnst_get_ind, cnst_mw
19 : use physconst, only: rga
20 : use oplus, only: oplus_init
21 : use edyn_init, only: edynamo_init
22 : use pio, only: var_desc_t
23 : use perf_mod, only: t_startf, t_stopf
24 : use epotential_params, only: epot_active, epot_crit_colats
25 : use shr_const_mod, only: SHR_CONST_REARTH ! meters
26 :
27 : implicit none
28 :
29 : private
30 :
31 : public :: ionosphere_readnl
32 : public :: ionosphere_init
33 : public :: ionosphere_run1
34 : public :: ionosphere_run2
35 : public :: ionosphere_init_restart
36 : public :: ionosphere_write_restart
37 : public :: ionosphere_read_restart
38 : public :: ionosphere_final
39 :
40 : ! private data
41 :
42 : ! opmmrtm1_phys is O+ at previous time step (phys grid decomposed)
43 : ! It needs to persist from time-step to time-step and across restarts
44 : ! On physics grid
45 : real(r8), allocatable :: opmmrtm1_phys(:,:,:)
46 : type(var_desc_t) :: Optm1_vdesc
47 : logical :: opmmrtm1_initialized
48 :
49 : integer :: index_ped, index_hall, index_te, index_ti
50 : integer :: index_ui, index_vi, index_wi
51 :
52 : integer :: ixo2=-1, ixo=-1, ixh=-1
53 : integer :: ixo2p=-1, ixnop=-1, ixn2p=-1, ixop=-1
54 :
55 : ! indices for accessing ions in pbuf when non-advected
56 : integer :: sIndxOp=-1, sIndxO2p=-1, sIndxNOp=-1, sIndxN2p=-1
57 :
58 : real(r8) :: rmassO2 ! O2 molecular weight kg/kmol
59 : real(r8) :: rmassO1 ! O atomic weight kg/kmol
60 : real(r8) :: rmassH ! H atomic weight kg/kmol
61 : real(r8) :: rmassN2 ! N2 molecular weight kg/kmol
62 : real(r8) :: rmassO2p ! O2+ molecular weight kg/kmol
63 : real(r8) :: rmassNOp ! NO+ molecular weight kg/kmol
64 : real(r8) :: rmassN2p ! N2+ molecular weight kg/kmol
65 : real(r8) :: rmassOp ! O+ molecular weight kg/kmol
66 :
67 : ! ionos_edyn_active == .true. will activate the edynamo which will
68 : ! generate ion drift velocities used in oplus transport, otherwise
69 : ! empirical ion drifts calculated in exbdrift (physics) will be used.
70 : logical, public, protected :: ionos_edyn_active = .true.
71 : logical, protected :: ionos_xport_active = .true. ! if true, call d_pie_coupling
72 : !
73 : logical, public, protected :: ionos_oplus_xport = .true. ! if true, call sub oplus (based on tiegcm oplus.F)
74 : integer, public, protected :: ionos_xport_nsplit = 5 ! number of substeps for O+ transport per model time step
75 : logical, public, protected :: oplus_ring_polar_filter = .false. ! switch to apply ring polar filter
76 :
77 : real(r8) :: oplus_adiff_limiter = 1.5e+8_r8 ! limiter for ambipolar diffusion coefficient
78 : real(r8) :: oplus_shapiro_const = 0.03_r8 ! shapiro constant for spatial smoother
79 : logical :: oplus_enforce_floor = .true. ! switch to apply Stan's floor
80 :
81 : integer, parameter :: max_num_files = 20
82 : character(len=cl) :: wei05_coefs_file = 'NONE' !'wei05sc.nc'
83 : character(len=cl) :: amienh_files(max_num_files) = 'NONE'
84 : character(len=cl) :: amiesh_files(max_num_files) = 'NONE'
85 : character(len=cl) :: ltr_files(max_num_files) = 'NONE'
86 :
87 :
88 : character(len=16) :: ionos_epotential_model = 'none'
89 : logical :: ionos_epotential_amie = .false.
90 : logical :: ionos_epotential_ltr = .false.
91 : integer :: indxefx=-1, indxkev=-1
92 :
93 : integer :: oplus_nlon, oplus_nlat ! Oplus grid
94 : integer :: ionos_npes = -1
95 :
96 : logical :: state_debug_checks = .false.
97 : logical :: ionos_debug_hist = .false.
98 :
99 : integer :: mag_nlon=0, mag_nlat=0, mag_nlev=0, mag_ngrid=0
100 :
101 : real(r8), parameter :: rearth_inv = 1._r8/SHR_CONST_REARTH ! /meters
102 :
103 : contains
104 :
105 : !---------------------------------------------------------------------------
106 : !---------------------------------------------------------------------------
107 768 : subroutine ionosphere_readnl( nlfile )
108 :
109 : use namelist_utils, only: find_group_name
110 : use units, only: getunit, freeunit
111 : use spmd_utils, only: masterproc, mpicom, masterprocid
112 : use spmd_utils, only: mpi_real8, mpi_logical, mpi_integer, mpi_character
113 : use cam_logfile, only: iulog
114 :
115 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
116 :
117 : ! Local variables
118 : integer :: unitn, ierr, ipos
119 : integer :: oplus_grid(2)
120 : character(len=8) :: edyn_grid
121 : integer :: total_pes
122 : character(len=*), parameter :: subname = 'ionosphere_readnl'
123 :
124 : namelist /ionosphere_nl/ ionos_xport_active, ionos_edyn_active, ionos_oplus_xport, ionos_xport_nsplit
125 : namelist /ionosphere_nl/ oplus_adiff_limiter, oplus_shapiro_const, oplus_enforce_floor, oplus_ring_polar_filter
126 : namelist /ionosphere_nl/ ionos_epotential_model, ionos_epotential_amie, ionos_epotential_ltr, wei05_coefs_file
127 : namelist /ionosphere_nl/ amienh_files, amiesh_files, wei05_coefs_file, ltr_files
128 : namelist /ionosphere_nl/ epot_crit_colats
129 : namelist /ionosphere_nl/ ionos_npes
130 : namelist /ionosphere_nl/ oplus_grid, edyn_grid
131 : namelist /ionosphere_nl/ ionos_debug_hist
132 :
133 768 : oplus_grid = 0
134 :
135 : ! Read namelist
136 768 : if (masterproc) then
137 2 : unitn = getunit()
138 2 : open( unitn, file=trim(nlfile), status='old' )
139 2 : call find_group_name(unitn, 'ionosphere_nl', status=ierr)
140 2 : if (ierr == 0) then
141 2 : read(unitn, ionosphere_nl, iostat=ierr)
142 2 : if (ierr /= 0) then
143 0 : call endrun(subname // ':: ERROR reading namelist')
144 : end if
145 : end if
146 2 : close(unitn)
147 2 : call freeunit(unitn)
148 : end if
149 :
150 : ! Broadcast namelist variables
151 768 : call mpi_bcast(ionos_xport_active, 1, mpi_logical, masterprocid, mpicom, ierr)
152 768 : call mpi_bcast(ionos_edyn_active, 1, mpi_logical, masterprocid, mpicom, ierr)
153 768 : call mpi_bcast(ionos_oplus_xport, 1, mpi_logical, masterprocid, mpicom, ierr)
154 768 : call mpi_bcast(ionos_xport_nsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
155 768 : call mpi_bcast(oplus_adiff_limiter, 1, mpi_real8, masterprocid, mpicom, ierr)
156 768 : call mpi_bcast(ionos_epotential_model, len(ionos_epotential_model), mpi_character, masterprocid, mpicom, ierr)
157 768 : call mpi_bcast(ionos_epotential_amie,1, mpi_logical, masterprocid, mpicom, ierr)
158 768 : call mpi_bcast(ionos_epotential_ltr,1, mpi_logical, masterprocid, mpicom, ierr)
159 768 : call mpi_bcast(wei05_coefs_file, len(wei05_coefs_file), mpi_character, masterprocid, mpicom, ierr)
160 768 : call mpi_bcast(amienh_files, max_num_files*len(amienh_files(1)), mpi_character, masterprocid, mpicom, ierr)
161 768 : call mpi_bcast(amiesh_files, max_num_files*len(amiesh_files(1)), mpi_character, masterprocid, mpicom, ierr)
162 768 : call mpi_bcast(ltr_files, max_num_files*len(ltr_files(1)), mpi_character, masterprocid, mpicom, ierr)
163 768 : call mpi_bcast(oplus_shapiro_const, 1, mpi_real8, masterprocid, mpicom, ierr)
164 768 : call mpi_bcast(oplus_enforce_floor, 1, mpi_logical, masterprocid, mpicom, ierr)
165 768 : call mpi_bcast(oplus_ring_polar_filter,1, mpi_logical, masterprocid, mpicom, ierr)
166 768 : call mpi_bcast(epot_crit_colats, 2, mpi_real8, masterprocid, mpicom, ierr)
167 768 : call mpi_bcast(ionos_npes, 1, mpi_integer, masterprocid, mpicom, ierr)
168 768 : call mpi_bcast(oplus_grid, 2, mpi_integer, masterprocid, mpicom, ierr)
169 768 : call mpi_bcast(edyn_grid, 8, mpi_character, masterprocid, mpicom, ierr)
170 768 : call mpi_bcast(ionos_debug_hist, 1, mpi_logical, masterprocid, mpicom, ierr)
171 :
172 : ! Extract grid settings
173 768 : oplus_nlon = oplus_grid(1)
174 768 : oplus_nlat = oplus_grid(2)
175 :
176 768 : ipos = scan(edyn_grid,'x')
177 768 : read(edyn_grid(:ipos-1),*) mag_nlon
178 768 : read(edyn_grid(ipos+1:),*) mag_nlat
179 :
180 768 : mag_nlev = 5 + int(log(real(mag_nlon,r8)/80._r8)/log(2._r8))
181 768 : mag_ngrid = (mag_nlon/10)*2
182 :
183 : ! Set npes in case of default settings
184 768 : call mpi_comm_size(mpicom, total_pes, ierr)
185 768 : if (ionos_npes<1) then
186 768 : ionos_npes = total_pes
187 0 : else if (ionos_npes>total_pes) then
188 0 : call endrun('ionosphere_readnl: ionos_npes > total_pes')
189 : end if
190 :
191 : ! log the user settings
192 768 : if (masterproc) then
193 2 : write(iulog,*) 'ionosphere_readnl: ionos_xport_active = ', ionos_xport_active
194 2 : write(iulog,*) 'ionosphere_readnl: ionos_edyn_active = ', ionos_edyn_active
195 2 : write(iulog,*) 'ionosphere_readnl: ionos_oplus_xport = ', ionos_oplus_xport
196 2 : write(iulog,*) 'ionosphere_readnl: ionos_xport_nsplit = ', ionos_xport_nsplit
197 2 : write(iulog,*) 'ionosphere_readnl: ionos_epotential_model = ', trim(ionos_epotential_model)
198 2 : write(iulog,*) 'ionosphere_readnl: ionos_epotential_amie = ', ionos_epotential_amie
199 2 : write(iulog,*) 'ionosphere_readnl: ionos_epotential_ltr = ', ionos_epotential_ltr
200 : write(iulog,'(a,2(g12.4))') &
201 2 : 'ionosphere_readnl: epot_crit_colats = ', epot_crit_colats
202 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: ionos_npes = ',ionos_npes
203 2 : write(iulog,*) 'ionosphere_readnl: oplus_adiff_limiter = ', oplus_adiff_limiter
204 2 : write(iulog,*) 'ionosphere_readnl: oplus_shapiro_const = ', oplus_shapiro_const
205 2 : write(iulog,*) 'ionosphere_readnl: oplus_enforce_floor = ', oplus_enforce_floor
206 2 : write(iulog,*) 'ionosphere_readnl: oplus_ring_polar_filter= ', oplus_ring_polar_filter
207 2 : if (ionos_xport_active) then
208 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: oplus_nlon = ',oplus_nlon
209 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: oplus_nlat = ',oplus_nlat
210 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: edyn_grid = '//edyn_grid
211 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: mag_nlon = ',mag_nlon
212 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: mag_nlat = ',mag_nlat
213 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: mag_nlev = ',mag_nlev
214 2 : write(iulog,'(a,i0)') 'ionosphere_readnl: mag_ngrid = ',mag_ngrid
215 : end if
216 : end if
217 768 : epot_active = .true.
218 :
219 768 : end subroutine ionosphere_readnl
220 :
221 : !---------------------------------------------------------------------------
222 : !---------------------------------------------------------------------------
223 768 : subroutine ionosphere_init()
224 : use spmd_utils, only: mpicom, iam
225 : use physics_buffer, only: pbuf_add_field, dtype_r8
226 : use cam_control_mod, only: initial_run
227 : use cam_history, only: addfld, add_default, horiz_only
228 : use edyn_mpi, only: mp_init
229 : use edyn_geogrid, only: set_geogrid
230 : use edyn_maggrid, only: alloc_maggrid
231 : use mo_apex, only: mo_apex_init1
232 : ! Hybrid level definitions:
233 : use ref_pres, only: pref_mid ! target alev(pver) midpoint levels
234 : use ref_pres, only: pref_edge ! target ailev(pverp) interface levels
235 : use amie_module, only: init_amie
236 : use ltr_module, only: init_ltr
237 : use wei05sc, only: weimer05_init
238 : use phys_control, only: phys_getopts
239 :
240 : ! local variables:
241 : integer :: sIndx
242 : character(len=*), parameter :: subname = 'ionosphere_init'
243 :
244 768 : call phys_getopts(state_debug_checks_out=state_debug_checks)
245 :
246 768 : if ( ionos_epotential_amie .or. ionos_epotential_ltr) then
247 0 : call pbuf_add_field('AUREFX', 'global', dtype_r8, (/pcols/), indxefx) ! Prescribed Energy flux
248 0 : call pbuf_add_field('AURKEV', 'global', dtype_r8, (/pcols/), indxkev) ! Prescribed Mean energy
249 : end if
250 768 : if (initial_run) then
251 : ! Read initial conditions (O+) on physics grid
252 384 : call ionosphere_read_ic()
253 : end if
254 :
255 768 : op_transport: if (ionos_xport_active) then
256 :
257 768 : index_ped = pbuf_get_index('PedConduct')
258 768 : index_hall = pbuf_get_index('HallConduct')
259 :
260 768 : index_te = pbuf_get_index('TElec')
261 768 : index_ti = pbuf_get_index('TIon')
262 : !
263 : ! pbuf indices to empirical ion drifts, to be passed to oplus_xport,
264 : ! if ionos_edyn_active is false.
265 : !
266 768 : index_ui = pbuf_get_index('UI')
267 768 : index_vi = pbuf_get_index('VI')
268 768 : index_wi = pbuf_get_index('WI')
269 :
270 : !---------------------------------------------------------------------
271 : ! Get indices for neutrals to get mixing ratios from state%q and masses
272 : !---------------------------------------------------------------------
273 768 : call cnst_get_ind('O2' ,ixo2 )
274 768 : call cnst_get_ind('O' ,ixo )
275 768 : call cnst_get_ind('H' ,ixh )
276 : !------------------------------------
277 : ! Get neutral molecular weights
278 : !------------------------------------
279 768 : rmassO2 = cnst_mw(ixo2)
280 768 : rmassO1 = cnst_mw(ixo)
281 768 : rmassH = cnst_mw(ixh)
282 768 : rmassN2 = 28._r8
283 :
284 768 : call cnst_get_ind('Op',ixop, abort=.false.)
285 768 : if (ixop > 0) then
286 0 : rMassOp = cnst_mw(ixop)
287 : else
288 768 : sIndxOp = slvd_index( 'Op' )
289 768 : if (sIndxOp > 0) then
290 768 : sIndx = get_spc_ndx( 'Op' )
291 768 : rmassOp = adv_mass(sIndx)
292 : else
293 0 : call endrun(subname//': Cannot find state or pbuf index for Op')
294 : end if
295 : end if
296 :
297 768 : call cnst_get_ind('O2p',ixo2p, abort=.false.)
298 768 : if (ixo2p > 0) then
299 0 : rMassO2p = cnst_mw(ixo2p)
300 : else
301 768 : sIndxO2p = slvd_index( 'O2p' )
302 768 : if (sIndxO2p > 0) then
303 768 : sIndx = get_spc_ndx( 'O2p' )
304 768 : rmassO2p = adv_mass(sIndx)
305 : else
306 0 : call endrun(subname//': Cannot find state or pbuf index for O2p')
307 : end if
308 : end if
309 :
310 768 : call cnst_get_ind('NOp',ixnop, abort=.false.)
311 768 : if (ixnop > 0) then
312 0 : rMassNOp = cnst_mw(ixnop)
313 : else
314 768 : sIndxNOp = slvd_index( 'NOp' )
315 768 : if (sIndxNOp > 0) then
316 768 : sIndx = get_spc_ndx( 'NOp' )
317 768 : rmassNOp = adv_mass(sIndx)
318 : else
319 0 : call endrun(subname//': Cannot find state or pbuf index for NOp')
320 : end if
321 : end if
322 :
323 768 : call cnst_get_ind('N2p',ixn2p, abort=.false.)
324 768 : if (ixn2p > 0) then
325 0 : rMassN2p = cnst_mw(ixn2p)
326 : else
327 768 : sIndxN2p = slvd_index( 'N2p' )
328 768 : if (sIndxN2p > 0) then
329 768 : sIndx = get_spc_ndx( 'N2p' )
330 768 : rmassN2p = adv_mass(sIndx)
331 : else
332 0 : call endrun(subname//': Cannot find state or pbuf index for N2p')
333 : end if
334 : end if
335 :
336 768 : call alloc_maggrid( mag_nlon, mag_nlat, mag_nlev, mag_ngrid )
337 :
338 768 : call mp_init(mpicom, ionos_npes, oplus_nlon, oplus_nlat, pver) ! set ntask,mytid
339 :
340 : ! set global geographic grid (sets coordinate distribution)
341 : ! lon0, lon1, etc. are set here
342 768 : call set_geogrid(oplus_nlon, oplus_nlat, pver, ionos_npes, iam, pref_mid, pref_edge)
343 :
344 768 : call edynamo_init(mpicom, ionos_debug_hist)
345 :
346 : call d_pie_init(ionos_edyn_active, ionos_oplus_xport, ionos_xport_nsplit, epot_crit_colats, &
347 768 : ionos_debug_hist)
348 :
349 768 : call ionosphere_alloc()
350 :
351 : call oplus_init(oplus_adiff_limiter, oplus_shapiro_const, oplus_enforce_floor, &
352 768 : oplus_ring_polar_filter, ionos_debug_hist)
353 :
354 1536 : call addfld('OpTM1&IC', (/ 'lev' /), 'I', 'kg/kg', 'O+ at time step minus 1', gridname='physgrid')
355 768 : call add_default ('OpTM1&IC',0, 'I')
356 :
357 : end if op_transport
358 :
359 : ! This has to be after edynamo_init (where maggrid is initialized)
360 768 : call mo_apex_init1()
361 :
362 768 : if (ionos_edyn_active) then
363 1536 : call addfld ('UI',(/ 'lev' /),'I','m/s', 'UI Zonal ion drift from edynamo')
364 1536 : call addfld ('VI',(/ 'lev' /),'I','m/s', 'VI Meridional ion drift from edynamo')
365 1536 : call addfld ('WI',(/ 'lev' /),'I','m/s', 'WI Vertical ion drift from edynamo')
366 1536 : call addfld ('UI&IC', (/ 'lev' /), 'I','m/s', 'Zonal ion drift velocity')
367 1536 : call addfld ('VI&IC', (/ 'lev' /), 'I','m/s', 'Meridional ion drift velocity')
368 1536 : call addfld ('WI&IC', (/ 'lev' /), 'I','m/s', 'Vertical ion drift velocity')
369 768 : call add_default ('UI&IC', 0, ' ')
370 768 : call add_default ('VI&IC', 0, ' ')
371 768 : call add_default ('WI&IC', 0, ' ')
372 : end if
373 768 : if ( ionos_epotential_amie ) then
374 0 : call init_amie(amienh_files,amiesh_files)
375 0 : call addfld ('amie_efx_phys', horiz_only, 'I', 'mW/m2', 'AMIE energy flux')
376 0 : call addfld ('amie_kev_phys', horiz_only, 'I', 'keV', 'AMIE mean energy')
377 : end if
378 768 : if ( ionos_epotential_ltr ) then
379 0 : call init_ltr(ltr_files)
380 0 : call addfld ('ltr_efx_phys', horiz_only, 'I', 'mW/m2', 'LTR energy flux')
381 0 : call addfld ('ltr_kev_phys', horiz_only, 'I', 'keV', 'LTR mean energy')
382 : end if
383 768 : if ( trim(ionos_epotential_model) == 'weimer' ) then
384 0 : call weimer05_init(wei05_coefs_file)
385 : end if
386 :
387 : ! d_pie_coupling diagnostics
388 : call addfld ('Z3GM', (/ 'lev' /), 'I', 'm', &
389 1536 : 'Geometric height', gridname='physgrid')
390 : call addfld ('Z3GMI', (/ 'lev' /), 'I', 'm', &
391 1536 : 'Geometric height (Interfaces)', gridname='physgrid')
392 :
393 768 : end subroutine ionosphere_init
394 :
395 : !----------------------------------------------------------------------------
396 : !----------------------------------------------------------------------------
397 8064 : subroutine ionosphere_run1(pbuf2d)
398 768 : use physics_buffer, only: physics_buffer_desc
399 : use cam_history, only: outfld, write_inithist
400 :
401 : ! args
402 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
403 :
404 : ! local vars
405 : integer :: i, j, lchnk, blksize ! indices
406 8064 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
407 :
408 8064 : real(r8), pointer :: pbuf_efx(:) ! Pointer to prescribed energy flux in pbuf
409 8064 : real(r8), pointer :: pbuf_kev(:) ! Pointer to prescribed mean energy in pbuf
410 :
411 : integer :: ncol
412 8064 : real(r8), pointer :: prescr_efx(:) ! prescribed energy flux
413 8064 : real(r8), pointer :: prescr_kev(:) ! prescribed characteristic mean energy
414 :
415 8064 : if( write_inithist() .and. ionos_xport_active ) then
416 0 : do lchnk = begchunk, endchunk
417 0 : call outfld ('OpTM1&IC', opmmrtm1_phys(:,:,lchnk), pcols, lchnk)
418 : end do
419 : end if
420 :
421 8064 : nullify(prescr_efx)
422 8064 : nullify(prescr_kev)
423 8064 : prescribed_epot: if ( ionos_epotential_amie .or. ionos_epotential_ltr ) then
424 0 : blksize = 0
425 0 : do lchnk = begchunk, endchunk
426 0 : blksize = blksize + get_ncols_p(lchnk)
427 : end do
428 :
429 0 : allocate(prescr_efx(blksize))
430 0 : allocate(prescr_kev(blksize))
431 :
432 : ! data assimilated potential
433 : call d_pie_epotent(ionos_epotential_model, epot_crit_colats, &
434 : cols=1, cole=blksize, efx_phys=prescr_efx, kev_phys=prescr_kev, &
435 0 : amie_in=ionos_epotential_amie, ltr_in=ionos_epotential_ltr )
436 :
437 : ! transform to pbuf for aurora...
438 :
439 0 : j = 0
440 0 : chnk_loop1: do lchnk = begchunk, endchunk
441 0 : ncol = get_ncols_p(lchnk)
442 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
443 0 : call pbuf_get_field(pbuf_chnk, indxefx, pbuf_efx)
444 0 : call pbuf_get_field(pbuf_chnk, indxkev, pbuf_kev)
445 :
446 0 : do i = 1, ncol
447 0 : j = j + 1
448 0 : pbuf_efx(i) = prescr_efx(j)
449 0 : pbuf_kev(i) = prescr_kev(j)
450 : end do
451 :
452 0 : if ( ionos_epotential_amie ) then
453 0 : call outfld('amie_efx_phys', pbuf_efx, pcols, lchnk)
454 0 : call outfld('amie_kev_phys', pbuf_kev, pcols, lchnk)
455 : endif
456 0 : if ( ionos_epotential_ltr) then
457 0 : call outfld('ltr_efx_phys', pbuf_efx, pcols, lchnk )
458 0 : call outfld('ltr_kev_phys', pbuf_kev, pcols, lchnk )
459 : end if
460 : end do chnk_loop1
461 :
462 0 : deallocate(prescr_efx, prescr_kev)
463 0 : nullify(prescr_efx)
464 : nullify(prescr_kev)
465 :
466 : else
467 :
468 : ! set cross tail potential before physics --
469 : ! aurora uses weimer derived potential
470 8064 : call d_pie_epotent( ionos_epotential_model, epot_crit_colats )
471 :
472 : end if prescribed_epot
473 :
474 16128 : end subroutine ionosphere_run1
475 :
476 : !---------------------------------------------------------------------------
477 : !---------------------------------------------------------------------------
478 14592 : subroutine ionosphere_run2(phys_state, pbuf2d)
479 :
480 8064 : use physics_types, only: physics_state
481 : use physics_buffer, only: physics_buffer_desc
482 : use cam_history, only: outfld, write_inithist, hist_fld_active
483 : use shr_assert_mod, only: shr_assert_in_domain
484 :
485 : ! - pull some fields from pbuf and dyn_in
486 : ! - invoke ionosphere/electro-dynamics coupling
487 : ! - push some fields back to physics via pbuf...
488 :
489 : ! args
490 : type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
491 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
492 :
493 : ! local vars
494 : integer :: i,j,k, lchnk
495 : integer :: astat
496 :
497 7296 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
498 :
499 7296 : real(r8), pointer :: sigma_ped_phys(:,:) ! Pedersen Conductivity from pbuf
500 7296 : real(r8), pointer :: sigma_hall_phys(:,:) ! Hall Conductivity from pbuf
501 7296 : real(r8), pointer :: te_phys(:,:) ! te from pbuf
502 7296 : real(r8), pointer :: ti_phys(:,:) ! ti from pbuf
503 7296 : real(r8), pointer :: mmrPO2p_phys(:,:) ! O2+ from pbuf
504 7296 : real(r8), pointer :: mmrPNOp_phys(:,:) ! NO+ from pbuf
505 7296 : real(r8), pointer :: mmrPN2p_phys(:,:) ! N2+ from pbuf
506 7296 : real(r8), pointer :: mmrPOp_phys(:,:) ! O+ from pbuf
507 : !
508 : ! Empirical ion drifts from exbdrift (to be converted to blocked for dpie_coupling):
509 7296 : real(r8), pointer :: ui_phys(:,:) ! zonal ion drift from pbuf
510 7296 : real(r8), pointer :: vi_phys(:,:) ! meridional ion drift from pbuf
511 7296 : real(r8), pointer :: wi_phys(:,:) ! vertical ion drift from pbuf
512 :
513 : integer :: ncol
514 :
515 : integer :: blksize ! number of columns in 2D block
516 :
517 7296 : real(r8), pointer :: sigma_ped_blck (:,:)
518 7296 : real(r8), pointer :: sigma_hall_blck(:,:)
519 7296 : real(r8), pointer :: ti_blck(:,:)
520 7296 : real(r8), pointer :: te_blck(:,:)
521 7296 : real(r8), pointer :: zi_blck(:,:) ! Geopotential on interfaces
522 7296 : real(r8), pointer :: hi_blck(:,:) ! Geometric height on interfaces
523 7296 : real(r8), pointer :: ui_blck(:,:)
524 7296 : real(r8), pointer :: vi_blck(:,:)
525 7296 : real(r8), pointer :: wi_blck(:,:)
526 7296 : real(r8), pointer :: omega_blck(:,:)
527 7296 : real(r8), pointer :: tn_blck(:,:)
528 :
529 : ! From physics state
530 7296 : real(r8), pointer :: u_blck(:,:)
531 7296 : real(r8), pointer :: v_blck(:,:)
532 7296 : real(r8), pointer :: pmid_blck(:,:)
533 7296 : real(r8), pointer :: phis(:) ! surface geopotential
534 : ! Constituents
535 7296 : real(r8), pointer :: n2mmr_blck(:,:)
536 7296 : real(r8), pointer :: o2mmr_blck(:,:)
537 7296 : real(r8), pointer :: o1mmr_blck(:,:)
538 7296 : real(r8), pointer :: h1mmr_blck(:,:)
539 7296 : real(r8), pointer :: o2pmmr_blck(:,:) ! O2+ (blocks)
540 7296 : real(r8), pointer :: nopmmr_blck(:,:) ! NO+ (blocks)
541 7296 : real(r8), pointer :: n2pmmr_blck(:,:) ! N2+ (blocks)
542 7296 : real(r8), pointer :: opmmr_blck(:,:) ! O+ (blocks)
543 7296 : real(r8), pointer :: opmmrtm1_blck(:,:) ! O+ previous time step (blocks)
544 7296 : real(r8), pointer :: mbar_blck(:,:) ! mean molecular weight
545 : ! Temp fields for outfld
546 : real(r8) :: r8tmp
547 : real(r8), pointer :: tempm(:,:) => null() ! Temp midpoint field for outfld
548 : real(r8), pointer :: tempi(:,:) => null() ! Temp interface field for outfld
549 : real(r8), parameter :: n2min = 1.e-6_r8 ! lower limit of N2 mixing ratios
550 :
551 : character(len=*), parameter :: subname = 'ionosphere_run2'
552 :
553 7296 : ionos_cpl: if (ionos_xport_active) then
554 :
555 7296 : blksize = 0
556 29184 : do lchnk = begchunk, endchunk
557 29184 : blksize = blksize + get_ncols_p(lchnk)
558 : end do
559 :
560 7296 : allocate(phis(pcols), stat=astat)
561 7296 : if (astat /= 0) then
562 0 : call endrun(subname//': failed to allocate phis')
563 : end if
564 21888 : allocate(u_blck(pver, blksize), stat=astat)
565 : if (astat /= 0) then
566 0 : call endrun(subname//': failed to allocate u_blck')
567 : end if
568 21888 : allocate(v_blck(pver, blksize), stat=astat)
569 : if (astat /= 0) then
570 0 : call endrun(subname//': failed to allocate v_blck')
571 : end if
572 21888 : allocate(sigma_ped_blck(pver, blksize), stat=astat)
573 : if (astat /= 0) then
574 0 : call endrun(subname//': failed to allocate sigma_ped_blck')
575 : end if
576 21888 : allocate(sigma_hall_blck(pver, blksize), stat=astat)
577 : if (astat /= 0) then
578 0 : call endrun(subname//': failed to allocate sigma_hall_blck')
579 : end if
580 21888 : allocate(ti_blck(pver, blksize), stat=astat)
581 : if (astat /= 0) then
582 0 : call endrun(subname//': failed to allocate ti_blck')
583 : end if
584 21888 : allocate(hi_blck(pver, blksize), stat=astat)
585 : if (astat /= 0) then
586 0 : call endrun(subname//': failed to allocate hi_blck')
587 : end if
588 21888 : allocate(te_blck(pver, blksize), stat=astat)
589 : if (astat /= 0) then
590 0 : call endrun(subname//': failed to allocate te_blck')
591 : end if
592 21888 : allocate(zi_blck(pver, blksize), stat=astat)
593 : if (astat /= 0) then
594 0 : call endrun(subname//': failed to allocate zi_blck')
595 : end if
596 21888 : allocate(ui_blck(pver, blksize), stat=astat)
597 : if (astat /= 0) then
598 0 : call endrun(subname//': failed to allocate ui_blck')
599 : end if
600 21888 : allocate(vi_blck(pver, blksize), stat=astat)
601 : if (astat /= 0) then
602 0 : call endrun(subname//': failed to allocate vi_blck')
603 : end if
604 21888 : allocate(wi_blck(pver, blksize), stat=astat)
605 : if (astat /= 0) then
606 0 : call endrun(subname//': failed to allocate wi_blck')
607 : end if
608 21888 : allocate(omega_blck(pver, blksize), stat=astat)
609 : if (astat /= 0) then
610 0 : call endrun(subname//': failed to allocate omega_blck')
611 : end if
612 21888 : allocate(tn_blck(pver, blksize), stat=astat)
613 : if (astat /= 0) then
614 0 : call endrun(subname//': failed to allocate tn_blck')
615 : end if
616 21888 : allocate(n2mmr_blck(pver, blksize), stat=astat)
617 : if (astat /= 0) then
618 0 : call endrun(subname//': failed to allocate n2mmr_blck')
619 : end if
620 21888 : allocate(o2mmr_blck(pver, blksize), stat=astat)
621 : if (astat /= 0) then
622 0 : call endrun(subname//': failed to allocate o2mmr_blck')
623 : end if
624 21888 : allocate(o1mmr_blck(pver, blksize), stat=astat)
625 : if (astat /= 0) then
626 0 : call endrun(subname//': failed to allocate o1mmr_blck')
627 : end if
628 21888 : allocate(h1mmr_blck(pver, blksize), stat=astat)
629 : if (astat /= 0) then
630 0 : call endrun(subname//': failed to allocate h1mmr_blck')
631 : end if
632 21888 : allocate(mbar_blck(pver, blksize), stat=astat)
633 : if (astat /= 0) then
634 0 : call endrun(subname//': failed to allocate mbar_blck')
635 : end if
636 21888 : allocate(pmid_blck(pver, blksize), stat=astat)
637 : if (astat /= 0) then
638 0 : call endrun(subname//': failed to allocate pmid_blck')
639 : end if
640 :
641 21888 : allocate(opmmrtm1_blck(pver, blksize), stat=astat)
642 : if (astat /= 0) then
643 0 : call endrun(subname//': failed to allocate opmmrtm1_blck')
644 : end if
645 :
646 7296 : if (sIndxOp > 0) then
647 21888 : allocate(opmmr_blck(pver, blksize), stat=astat)
648 : if (astat /= 0) then
649 0 : call endrun(subname//': failed to allocate opmmr_blck')
650 : end if
651 : end if
652 7296 : if (sIndxO2p > 0) then
653 21888 : allocate(o2pmmr_blck(pver, blksize), stat=astat)
654 : if (astat /= 0) then
655 0 : call endrun(subname//': failed to allocate o2pmmr_blck')
656 : end if
657 : end if
658 7296 : if (sIndxNOp > 0) then
659 21888 : allocate(nopmmr_blck(pver, blksize), stat=astat)
660 : if (astat /= 0) then
661 0 : call endrun(subname//': failed to allocate nopmmr_blck')
662 : end if
663 : end if
664 7296 : if (sIndxN2p > 0) then
665 21888 : allocate(n2pmmr_blck(pver, blksize), stat=astat)
666 : if (astat /= 0) then
667 0 : call endrun(subname//': failed to allocate n2pmmr_blck')
668 : end if
669 : end if
670 :
671 7296 : if (hist_fld_active('Z3GM')) then
672 7296 : allocate(tempm(pcols, pver))
673 : end if
674 :
675 7296 : if (hist_fld_active('Z3GMI')) then
676 0 : allocate(tempi(pcols, pver))
677 : end if
678 :
679 7296 : if (.not.opmmrtm1_initialized) then
680 0 : do lchnk = begchunk, endchunk
681 0 : ncol = get_ncols_p(lchnk)
682 :
683 0 : if (sIndxOp > 0) then
684 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
685 0 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPOp_phys, start=(/1,1,sIndxOp/), kount=(/pcols,pver,1/) )
686 0 : opmmrtm1_phys(:ncol,:pver,lchnk) = mmrPOp_phys(:ncol,:pver)
687 : else
688 0 : opmmrtm1_phys(:ncol,:pver,lchnk) = phys_state(lchnk)%q(:ncol,:pver, ixop)
689 : endif
690 : enddo
691 0 : opmmrtm1_initialized=.true.
692 : endif
693 :
694 7296 : j = 0
695 29184 : do lchnk = begchunk, endchunk
696 21888 : ncol = get_ncols_p(lchnk)
697 21888 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
698 :
699 : ! Gather data stored in pbuf and collect into blocked arrays
700 : ! Get Pedersen and Hall conductivities:
701 21888 : call pbuf_get_field(pbuf_chnk, index_ped, sigma_ped_phys)
702 21888 : call pbuf_get_field(pbuf_chnk, index_hall, sigma_hall_phys)
703 : ! Get ion and electron temperatures
704 21888 : call pbuf_get_field(pbuf_chnk, index_te, te_phys)
705 21888 : call pbuf_get_field(pbuf_chnk, index_ti, ti_phys)
706 : ! Get components of ion drift velocities
707 21888 : call pbuf_get_field(pbuf_chnk, index_ui, ui_phys)
708 21888 : call pbuf_get_field(pbuf_chnk, index_vi, vi_phys)
709 21888 : call pbuf_get_field(pbuf_chnk, index_wi, wi_phys)
710 : !--------------------------------------------------------
711 : ! Get ions from physics buffer if non-transported
712 : !--------------------------------------------------------
713 21888 : if (sIndxO2p > 0) then
714 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPO2p_phys, &
715 87552 : start=(/1,1,sIndxO2p/), kount=(/pcols,pver,1/) )
716 : end if
717 21888 : if (sIndxNOp > 0) then
718 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPNOp_phys, &
719 87552 : start=(/1,1,sIndxNOp/), kount=(/pcols,pver,1/) )
720 : end if
721 21888 : if (sIndxN2p > 0) then
722 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPN2p_phys, &
723 87552 : start=(/1,1,sIndxN2p/), kount=(/pcols,pver,1/) )
724 : end if
725 21888 : if (sIndxOp > 0) then
726 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPOp_phys, &
727 87552 : start=(/1,1,sIndxOp/), kount=(/pcols,pver,1/) )
728 : end if
729 :
730 : ! PHIS is from physics state
731 284544 : phis(:ncol) = phys_state(lchnk)%phis(:ncol)
732 284544 : do i = 1, ncol
733 262656 : j = j + 1
734 34429824 : do k = 1, pver
735 : ! physics state fields on levels
736 34145280 : u_blck(k, j) = phys_state(lchnk)%u(i, k)
737 34145280 : v_blck(k, j) = phys_state(lchnk)%v(i, k)
738 : !------------------------------------------------------------
739 : ! Might need geometric height on midpoints for output
740 : !------------------------------------------------------------
741 34145280 : if (hist_fld_active('Z3GM')) then
742 : ! geometric altitude (meters above sea level)
743 34145280 : tempm(i,k) = geometric_hgt(zgp=phys_state(lchnk)%zm(i,k), zsf=phis(i)*rga)
744 : end if
745 : ! physics state fields on interfaces (but only to pver)
746 34145280 : zi_blck(k, j) = phys_state(lchnk)%zi(i, k) + phis(i)*rga
747 : !------------------------------------------------------------
748 : ! Convert geopotential to geometric height at interfaces:
749 : !------------------------------------------------------------
750 : ! Note: zht is pver instead of pverp because dynamo does not
751 : ! use bottom interface
752 34145280 : hi_blck(k,j) = geometric_hgt(zgp=phys_state(lchnk)%zi(i,k), zsf=phis(i)*rga)
753 34145280 : if (hist_fld_active('Z3GMI')) then
754 0 : tempi(i,k) = hi_blck(k, j)
755 : end if
756 34145280 : omega_blck(k, j) = phys_state(lchnk)%omega(i, k)
757 34145280 : tn_blck(k, j) = phys_state(lchnk)%t(i, k)
758 34145280 : pmid_blck(k, j) = phys_state(lchnk)%pmid(i, k)
759 : ! Pedersen and Hall conductivities:
760 34145280 : sigma_ped_blck(k, j) = sigma_ped_phys(i, k)
761 34145280 : sigma_hall_blck(k, j) = sigma_hall_phys(i, k)
762 : ! ion and electron temperatures
763 34145280 : te_blck(k, j) = te_phys(i, k)
764 34145280 : ti_blck(k, j) = ti_phys(i, k)
765 : ! components of ion drift velocities
766 34145280 : ui_blck(k, j) = ui_phys(i, k)
767 34145280 : vi_blck(k, j) = vi_phys(i, k)
768 34145280 : wi_blck(k, j) = wi_phys(i, k)
769 : !------------------------------------------------------------
770 : ! ions from physics state if transported, otherwise from pbuf
771 : !------------------------------------------------------------
772 34145280 : if (ixo2p > 0) then
773 0 : o2pmmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixo2p)
774 34145280 : else if (sIndxO2p > 0) then
775 34145280 : o2pmmr_blck(k, j) = mmrPO2p_phys(i, k)
776 : else
777 0 : call endrun(subname//': No source for O2p')
778 : end if
779 34145280 : if (ixnop > 0) then
780 0 : nopmmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixnop)
781 34145280 : else if (sIndxNOp > 0) then
782 34145280 : nopmmr_blck(k, j) = mmrPNOp_phys(i, k)
783 : else
784 0 : call endrun(subname//': No source for NOp')
785 : end if
786 34145280 : if (ixn2p > 0) then
787 0 : n2pmmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixn2p)
788 34145280 : else if (sIndxN2p > 0) then
789 34145280 : n2pmmr_blck(k, j) = mmrPN2p_phys(i, k)
790 : else
791 0 : call endrun(subname//': No source for N2p')
792 : end if
793 34145280 : if (ixop > 0) then
794 0 : opmmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixop)
795 34145280 : else if (sIndxOp > 0) then
796 34145280 : opmmr_blck(k, j) = mmrPOp_phys(i, k)
797 : else
798 0 : call endrun(subname//': No source for Op')
799 : end if
800 34145280 : opmmrtm1_blck(k, j) = opmmrtm1_phys(i, k, lchnk)
801 : !------------------------------------
802 : ! neutrals from advected tracers array
803 : !------------------------------------
804 34145280 : o2mmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixo2)
805 34145280 : o1mmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixo)
806 34407936 : h1mmr_blck(k, j) = phys_state(lchnk)%q(i, k, ixh)
807 : end do
808 : end do ! do i = 1, ncol
809 :
810 : !------------------------------------------------------------------
811 : ! Save OMEGA and analytically derived geometric height
812 : !------------------------------------------------------------------
813 21888 : if (hist_fld_active('Z3GM')) then
814 14249088 : tempm(ncol+1:, :) = 0.0_r8
815 21888 : call outfld('Z3GM', tempm, pcols, lchnk)
816 : end if
817 29184 : if (hist_fld_active('Z3GMI')) then
818 0 : tempi(ncol+1:, :) = 0.0_r8
819 0 : call outfld('Z3GMI', tempi, pcols, lchnk)
820 : end if
821 : end do ! do lchnk = begchunk, endchunk
822 :
823 : !---------------------------------------------------------------------
824 : ! Compute and save mean molecular weight:
825 : !---------------------------------------------------------------------
826 7296 : j = 0
827 29184 : do lchnk = begchunk, endchunk
828 21888 : ncol = get_ncols_p(lchnk)
829 291840 : do i = 1, ncol
830 262656 : j = j + 1
831 34429824 : do k = 1, pver
832 34145280 : r8tmp = o1mmr_blck(k,j) + o2mmr_blck(k,j) + h1mmr_blck(k,j)
833 34145280 : n2mmr_blck(k, j) = max(1.0_r8 - r8tmp, n2min)
834 34145280 : r8tmp = o1mmr_blck(k, j) / rmassO1
835 34145280 : r8tmp = r8tmp + (o2mmr_blck(k, j) / rmassO2)
836 34145280 : r8tmp = r8tmp + (h1mmr_blck(k, j) / rmassH)
837 34145280 : r8tmp = r8tmp + (n2mmr_blck(k, j) / rmassN2)
838 34407936 : mbar_blck(k, j) = 1.0_r8 / r8tmp
839 : end do
840 : end do
841 : end do
842 :
843 7296 : call t_startf('d_pie_coupling')
844 :
845 : ! Compute geometric height and some diagnostic fields needed by
846 : ! the dynamo. Output some fields from physics grid
847 : ! This code is inside the timer as it is part of the coupling
848 : !
849 : ! waccmx ionosphere electro-dynamics -- transports O+ and
850 : ! provides updates to ion drift velocities (on physics grid)
851 : ! All fields are on physics mesh, (pver, blksize),
852 : ! where blksize is the total number of columns on this task
853 :
854 : call d_pie_coupling(omega_blck, pmid_blck, zi_blck, hi_blck, &
855 : u_blck, v_blck, tn_blck, sigma_ped_blck, sigma_hall_blck, &
856 : te_blck, ti_blck, mbar_blck, n2mmr_blck, o2mmr_blck, &
857 : o1mmr_blck, o2pmmr_blck, nopmmr_blck, n2pmmr_blck, &
858 : opmmr_blck, opmmrtm1_blck, ui_blck, vi_blck, wi_blck, &
859 7296 : rmassO2p, rmassNOp, rmassN2p, rmassOp, 1, blksize, pver)
860 :
861 7296 : call t_stopf ('d_pie_coupling')
862 :
863 7296 : if (state_debug_checks) then
864 7296 : call shr_assert_in_domain(ui_blck, is_nan=.false., varname="ui_blck", msg="NaN found in ionosphere_run2")
865 7296 : call shr_assert_in_domain(vi_blck, is_nan=.false., varname="vi_blck", msg="NaN found in ionosphere_run2")
866 7296 : call shr_assert_in_domain(wi_blck, is_nan=.false., varname="wi_blck", msg="NaN found in ionosphere_run2")
867 7296 : call shr_assert_in_domain(opmmr_blck, is_nan=.false., varname="opmmr_blck", msg="NaN found in ionosphere_run2")
868 : end if
869 :
870 : !
871 : !----------------------------------------
872 : ! Put data back in to state or pbuf
873 : !----------------------------------------
874 : ! blocks --> physics chunks
875 :
876 7296 : j = 0
877 29184 : do lchnk = begchunk, endchunk
878 21888 : ncol = phys_state(lchnk)%ncol
879 21888 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
880 :
881 21888 : call pbuf_get_field(pbuf_chnk, index_ui, ui_phys)
882 21888 : call pbuf_get_field(pbuf_chnk, index_vi, vi_phys)
883 21888 : call pbuf_get_field(pbuf_chnk, index_wi, wi_phys)
884 21888 : if (sIndxOp > 0) then
885 : call pbuf_get_field(pbuf_chnk, slvd_pbf_ndx, mmrPOp_phys, &
886 87552 : start=(/1,1,sIndxOp/), kount=(/pcols,pver,1/))
887 : end if
888 284544 : do i = 1, ncol
889 262656 : j = j + 1
890 34429824 : do k = 1, pver
891 34145280 : ui_phys(i, k) = ui_blck(k, j)
892 34145280 : vi_phys(i, k) = vi_blck(k, j)
893 34145280 : wi_phys(i, k) = wi_blck(k, j)
894 34145280 : if (ixop > 0) then
895 0 : phys_state(lchnk)%q(i, k, ixop) = opmmr_blck(k, j)
896 34145280 : else if (sIndxOp > 0) then
897 34145280 : mmrPOp_phys(i, k) = opmmr_blck(k, j)
898 : else
899 0 : call endrun(subname//': No destination for Op')
900 : end if
901 34407936 : opmmrtm1_phys(i,k,lchnk) = opmmrtm1_blck(k,j)
902 : end do
903 : end do
904 :
905 29184 : if (ionos_edyn_active) then
906 21888 : call outfld('UI', ui_phys, pcols, lchnk)
907 21888 : call outfld('VI', vi_phys, pcols, lchnk)
908 21888 : call outfld('WI', wi_phys, pcols, lchnk)
909 21888 : if (write_inithist()) then
910 0 : call outfld('UI&IC', ui_phys, pcols, lchnk)
911 0 : call outfld('VI&IC', vi_phys, pcols, lchnk)
912 0 : call outfld('WI&IC', wi_phys, pcols, lchnk)
913 : end if
914 : end if
915 : end do
916 :
917 7296 : if (associated(opmmr_blck)) then
918 7296 : deallocate(opmmr_blck)
919 : nullify(opmmr_blck)
920 : end if
921 7296 : if (associated(o2pmmr_blck)) then
922 7296 : deallocate(o2pmmr_blck)
923 : nullify(o2pmmr_blck)
924 : end if
925 7296 : if (associated(nopmmr_blck)) then
926 7296 : deallocate(nopmmr_blck)
927 : nullify(nopmmr_blck)
928 : end if
929 7296 : if (associated(n2pmmr_blck)) then
930 7296 : deallocate(n2pmmr_blck)
931 : nullify(n2pmmr_blck)
932 : end if
933 7296 : if (associated(tempi)) then
934 0 : deallocate(tempi)
935 : nullify(tempi)
936 : end if
937 7296 : if (associated(tempm)) then
938 7296 : deallocate(tempm)
939 : nullify(tempm)
940 : end if
941 7296 : deallocate(opmmrtm1_blck)
942 : nullify(opmmrtm1_blck)
943 7296 : deallocate(phis)
944 : nullify(phis)
945 7296 : deallocate(u_blck)
946 : nullify(u_blck)
947 7296 : deallocate(v_blck)
948 : nullify(v_blck)
949 7296 : deallocate(sigma_ped_blck)
950 : nullify(sigma_ped_blck)
951 7296 : deallocate(sigma_hall_blck)
952 : nullify(sigma_hall_blck)
953 7296 : deallocate(ti_blck)
954 : nullify(ti_blck)
955 7296 : deallocate(hi_blck)
956 : nullify(hi_blck)
957 7296 : deallocate(te_blck)
958 : nullify(te_blck)
959 7296 : deallocate(zi_blck)
960 : nullify(zi_blck)
961 7296 : deallocate(ui_blck)
962 : nullify(ui_blck)
963 7296 : deallocate(vi_blck)
964 : nullify(vi_blck)
965 7296 : deallocate(wi_blck)
966 : nullify(wi_blck)
967 7296 : deallocate(omega_blck)
968 : nullify(omega_blck)
969 7296 : deallocate(tn_blck)
970 : nullify(tn_blck)
971 7296 : deallocate(n2mmr_blck)
972 : nullify(n2mmr_blck)
973 7296 : deallocate(o2mmr_blck)
974 : nullify(o2mmr_blck)
975 7296 : deallocate(o1mmr_blck)
976 : nullify(o1mmr_blck)
977 7296 : deallocate(h1mmr_blck)
978 : nullify(h1mmr_blck)
979 7296 : deallocate(mbar_blck)
980 : nullify(mbar_blck)
981 7296 : deallocate(pmid_blck)
982 : nullify(pmid_blck)
983 :
984 : end if ionos_cpl
985 :
986 14592 : end subroutine ionosphere_run2
987 :
988 : !---------------------------------------------------------------------------
989 : !---------------------------------------------------------------------------
990 768 : subroutine ionosphere_init_restart(File)
991 7296 : use pio, only: file_desc_t, pio_double, pio_def_var
992 : use cam_pio_utils, only: cam_pio_def_dim
993 : use cam_grid_support, only: cam_grid_id, cam_grid_write_attr
994 : use cam_grid_support, only: cam_grid_header_info_t
995 :
996 : type(File_desc_t), intent(inout) :: File
997 :
998 : integer :: grid_id
999 : integer :: hdimcnt, ierr, i
1000 : integer :: dimids(3), ndims
1001 768 : type(cam_grid_header_info_t) :: info
1002 :
1003 768 : if (ionos_xport_active) then
1004 768 : grid_id = cam_grid_id('physgrid')
1005 768 : call cam_grid_write_attr(File, grid_id, info)
1006 768 : hdimcnt = info%num_hdims()
1007 2304 : do i = 1, hdimcnt
1008 2304 : dimids(i) = info%get_hdimid(i)
1009 : end do
1010 768 : ndims = hdimcnt + 1
1011 :
1012 768 : call cam_pio_def_dim(File, 'lev', pver, dimids(ndims), &
1013 1536 : existOK=.true.)
1014 :
1015 : ierr = pio_def_var(File, 'Optm1', pio_double, dimids(1:ndims), &
1016 768 : Optm1_vdesc)
1017 : end if
1018 768 : end subroutine ionosphere_init_restart
1019 :
1020 : !---------------------------------------------------------------------------
1021 : !---------------------------------------------------------------------------
1022 768 : subroutine ionosphere_write_restart(File)
1023 768 : use pio, only: io_desc_t, file_desc_t, pio_write_darray
1024 : use pio, only: pio_double
1025 : use cam_grid_support, only: cam_grid_id, cam_grid_write_var
1026 : use cam_grid_support, only: cam_grid_get_decomp, cam_grid_dimensions
1027 : use phys_grid, only: phys_decomp
1028 :
1029 : type(file_desc_t), intent(inout) :: File
1030 :
1031 : integer :: ierr
1032 : integer :: physgrid
1033 : integer :: dims(3), gdims(3)
1034 : integer :: nhdims
1035 : type(io_desc_t), pointer :: iodesc3d
1036 :
1037 768 : if (ionos_xport_active) then
1038 :
1039 : ! Write grid vars
1040 768 : call cam_grid_write_var(File, phys_decomp)
1041 :
1042 768 : physgrid = cam_grid_id('physgrid')
1043 768 : call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)
1044 768 : nhdims = nhdims + 1
1045 768 : gdims(nhdims) = pver
1046 768 : dims(1) = pcols
1047 768 : dims(2) = pver
1048 768 : dims(3) = endchunk - begchunk + 1
1049 : call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), &
1050 768 : pio_double, iodesc3d)
1051 :
1052 768 : call pio_write_darray(File, Optm1_vdesc, iodesc3d, opmmrtm1_phys, ierr)
1053 : end if
1054 :
1055 768 : end subroutine ionosphere_write_restart
1056 :
1057 : !---------------------------------------------------------------------------
1058 : !---------------------------------------------------------------------------
1059 384 : subroutine ionosphere_read_restart(File)
1060 768 : use pio, only: io_desc_t, file_desc_t, pio_inq_varid
1061 : use pio, only: pio_read_darray, pio_double
1062 : use cam_grid_support, only: cam_grid_id
1063 : use cam_grid_support, only: cam_grid_get_decomp, cam_grid_dimensions
1064 :
1065 : type(file_desc_t), intent(inout) :: File
1066 :
1067 : integer :: ierr
1068 : integer :: physgrid
1069 : integer :: dims(3), gdims(3)
1070 : integer :: nhdims
1071 : type(io_desc_t), pointer :: iodesc3d
1072 :
1073 384 : if (ionos_xport_active) then
1074 384 : call ionosphere_alloc()
1075 :
1076 384 : physgrid = cam_grid_id('physgrid')
1077 384 : call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)
1078 384 : nhdims = nhdims + 1
1079 384 : gdims(nhdims) = pver
1080 384 : dims(1) = pcols
1081 384 : dims(2) = pver
1082 384 : dims(3) = endchunk - begchunk + 1
1083 : call cam_grid_get_decomp(physgrid, dims(1:3), gdims(1:nhdims), &
1084 384 : pio_double, iodesc3d)
1085 :
1086 384 : ierr = pio_inq_varid(File, 'Optm1', Optm1_vdesc)
1087 384 : call pio_read_darray(File, Optm1_vdesc, iodesc3d, opmmrtm1_phys, ierr)
1088 384 : opmmrtm1_initialized = .true.
1089 : end if
1090 :
1091 384 : end subroutine ionosphere_read_restart
1092 :
1093 : !---------------------------------------------------------------------------
1094 : !---------------------------------------------------------------------------
1095 768 : subroutine ionosphere_final
1096 :
1097 384 : use edyn_esmf, only: edyn_esmf_final
1098 :
1099 768 : call edyn_esmf_final()
1100 :
1101 768 : if (allocated(opmmrtm1_phys)) then
1102 768 : deallocate(opmmrtm1_phys)
1103 : end if
1104 :
1105 768 : end subroutine ionosphere_final
1106 :
1107 : !===========================================================================
1108 : !---------------------------------------------------------------------------
1109 : !---------------------------------------------------------------------------
1110 384 : subroutine ionosphere_read_ic()
1111 :
1112 768 : use pio, only: file_desc_t
1113 : use ncdio_atm, only: infld
1114 : use cam_initfiles, only: initial_file_get_id
1115 : use cam_grid_support, only: cam_grid_check, cam_grid_id
1116 : use cam_grid_support, only: cam_grid_get_dim_names
1117 :
1118 : type(file_desc_t), pointer :: fh_ini ! PIO filehandle
1119 :
1120 : integer :: grid_id ! grid ID for data mapping
1121 : character(len=8) :: dim1name, dim2name
1122 : logical :: readvar
1123 : character(len=*), parameter :: subname = 'ionosphere_read_ic'
1124 :
1125 384 : if ( ionos_xport_active ) then
1126 384 : call ionosphere_alloc()
1127 :
1128 384 : fh_ini => initial_file_get_id()
1129 384 : grid_id = cam_grid_id('physgrid')
1130 384 : if (.not. cam_grid_check(grid_id)) then
1131 0 : call endrun(trim(subname)//': Internal error, no "physgrid" grid')
1132 : end if
1133 384 : call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
1134 :
1135 : ! try reading in OpTM1 from the IC file
1136 : call infld('OpTM1', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, &
1137 384 : begchunk, endchunk, opmmrtm1_phys, readvar, gridname='physgrid')
1138 384 : if (.not. readvar) then
1139 : ! if OpTM1 is not included in the IC file then try using O+
1140 : call infld('Op', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, &
1141 0 : begchunk, endchunk, opmmrtm1_phys, readvar, gridname='physgrid')
1142 : end if
1143 384 : opmmrtm1_initialized = readvar
1144 : end if
1145 :
1146 384 : end subroutine ionosphere_read_ic
1147 :
1148 : !---------------------------------------------------------------------------
1149 : !---------------------------------------------------------------------------
1150 1536 : subroutine ionosphere_alloc()
1151 384 : use infnan, only: nan, assignment(=)
1152 : integer :: astat
1153 :
1154 1536 : if (.not. allocated(opmmrtm1_phys)) then
1155 2304 : allocate(opmmrtm1_phys(pcols, pver, begchunk:endchunk), stat=astat)
1156 768 : if (astat /= 0) then
1157 0 : call endrun('ionosphere_alloc: failed to allocate opmmrtm1_phys')
1158 : end if
1159 768 : opmmrtm1_phys = nan
1160 768 : opmmrtm1_initialized = .false.
1161 : end if
1162 :
1163 1536 : end subroutine ionosphere_alloc
1164 :
1165 : !==========================================================================
1166 :
1167 : ! calculates geometric height (meters above sea level)
1168 68290560 : pure function geometric_hgt( zgp, zsf ) result(zgm)
1169 :
1170 : real(r8), intent(in) :: zgp ! geopotential height (m)
1171 : real(r8), intent(in) :: zsf ! surface height above sea level (m)
1172 : real(r8) :: zgm ! geometric height above sea level (m)
1173 :
1174 : real(r8) :: tmp
1175 :
1176 : ! Hanli's formulation:
1177 : ! Z_gm = 1/(1 - (1+Zs/r) * Z_gp/r) * (Zs + (1+Zs/r) * Z_gp)
1178 : ! Z_gm: geometric height
1179 : ! Zs: Surface height
1180 : ! Z_gp: model calculated geopotential height (zm and zi in the model)
1181 :
1182 68290560 : tmp = 1._r8+zsf*rearth_inv
1183 68290560 : zgm = (zsf + tmp*zgp) / (1._r8 - tmp*zgp*rearth_inv)
1184 :
1185 68290560 : end function geometric_hgt
1186 :
1187 : end module ionosphere_interface
|
