Line data Source code
1 : module dpie_coupling
2 : !
3 : ! Dynamics/Physics Ionosphere/Electrodynamics coupler.
4 : !
5 : use shr_kind_mod, only: r8 => shr_kind_r8
6 : use cam_logfile, only: iulog
7 : use cam_history, only: hist_fld_active, outfld
8 : use cam_history, only: addfld, horiz_only
9 : use cam_history_support, only: fillvalue
10 : use cam_abortutils, only: endrun
11 : use spmd_utils, only: masterproc, mpi_logical, mpicom, masterprocid
12 : use edyn_mpi, only: array_ptr_type
13 : use perf_mod, only: t_startf, t_stopf
14 : use amie_module, only: getamie
15 : use ltr_module, only: getltr
16 : use edyn_solve, only: phihm
17 : use edyn_params, only: dtr, rtd
18 : use aurora_params, only: prescribed_period ! turns on overwrite of energy fields in aurora phys
19 :
20 : implicit none
21 :
22 : private
23 : public :: d_pie_init
24 : public :: d_pie_epotent ! sets electric potential
25 : public :: d_pie_coupling ! handles coupling with edynamo and ion transport
26 :
27 : logical :: ionos_edyn_active ! if true, call oplus_xport for O+ transport
28 : logical :: ionos_oplus_xport ! if true, call oplus_xport for O+ transport
29 : integer :: nspltop ! nsplit for oplus_xport
30 :
31 : logical :: debug = .false.
32 :
33 : real(r8) :: crad(2), crit(2)
34 : logical :: crit_user_set = .false.
35 : real(r8), parameter :: amie_default_crit(2) = (/ 35._r8, 40._r8 /)
36 :
37 : logical :: debug_hist
38 :
39 : contains
40 : !----------------------------------------------------------------------
41 768 : subroutine d_pie_init( edyn_active_in, oplus_xport_in, oplus_nsplit_in, crit_colats_deg, ionos_debug_hist )
42 :
43 : logical, intent(in) :: edyn_active_in, oplus_xport_in
44 : integer, intent(in) :: oplus_nsplit_in
45 : real(r8),intent(in) :: crit_colats_deg(:)
46 : logical, intent(in) :: ionos_debug_hist
47 :
48 768 : debug_hist = ionos_debug_hist
49 :
50 768 : ionos_edyn_active = edyn_active_in
51 768 : ionos_oplus_xport = oplus_xport_in
52 768 : nspltop = oplus_nsplit_in
53 :
54 768 : crit_user_set = all( crit_colats_deg(:) > 0._r8 )
55 768 : if (crit_user_set) then
56 0 : crit(:) = crit_colats_deg(:)*dtr
57 : end if
58 :
59 768 : call addfld ('HMF2' , horiz_only , 'I', 'km' ,'Height of the F2 Layer' , gridname='physgrid')
60 768 : call addfld ('NMF2' , horiz_only , 'I', 'cm-3','Peak Density of the F2 Layer', gridname='physgrid')
61 :
62 1536 : call addfld ('OpDens' ,(/ 'lev' /), 'I', 'cm^-3','O+ Number Density' , gridname='physgrid')
63 1536 : call addfld ('EDens' ,(/ 'lev' /), 'I', 'cm^-3','e Number Density (sum of O2+,NO+,N2+,O+)', gridname='physgrid')
64 :
65 768 : call addfld ('prescr_efxp' , horiz_only, 'I','mW/m2','Prescribed energy flux on geo grid' ,gridname='physgrid')
66 768 : call addfld ('prescr_kevp' , horiz_only, 'I','keV ','Prescribed mean energy on geo grid' ,gridname='physgrid')
67 :
68 768 : call addfld ('prescr_phihm' , horiz_only, 'I','VOLTS','Prescribed Electric Potential-mag grid' ,gridname='gmag_grid')
69 768 : call addfld ('prescr_efxm' , horiz_only, 'I','mW/m2','Prescribed energy flux on mag grid' ,gridname='gmag_grid')
70 768 : call addfld ('prescr_kevm' , horiz_only, 'I','keV ','Prescribed mean energy on mag grid' ,gridname='gmag_grid')
71 :
72 768 : if (debug_hist) then
73 : ! Dynamo inputs (called from dpie_coupling. Fields are in waccm format, in CGS units):
74 0 : call addfld ('DPIE_OMEGA',(/ 'lev' /), 'I', 'Pa/s ','OMEGA input to DPIE coupling', gridname='physgrid')
75 0 : call addfld ('DPIE_MBAR' ,(/ 'lev' /), 'I', 'kg/kmole','MBAR Mean Mass from dpie_coupling', gridname='physgrid')
76 0 : call addfld ('DPIE_TN ',(/ 'lev' /), 'I', 'deg K ','DPIE_TN' , gridname='physgrid')
77 0 : call addfld ('DPIE_UN ',(/ 'lev' /), 'I', 'cm/s ','DPIE_UN' , gridname='physgrid')
78 0 : call addfld ('DPIE_VN ',(/ 'lev' /), 'I', 'cm/s ','DPIE_VN' , gridname='physgrid')
79 0 : call addfld ('DPIE_WN ',(/ 'lev' /), 'I', 'cm/s ','DPIE_WN' , gridname='physgrid')
80 0 : call addfld ('DPIE_OM ',(/ 'lev' /), 'I', 's-1 ','DPIE_OM' , gridname='physgrid')
81 0 : call addfld ('DPIE_ZHT ',(/ 'lev' /), 'I', 'cm ','DPIE_ZHT (geometric height,simple)', gridname='physgrid')
82 0 : call addfld ('DPIE_ZGI ',(/ 'lev' /), 'I', 'cm ','DPIE_ZGI (geopotential height on interfaces)', gridname='physgrid')
83 0 : call addfld ('DPIE_O2 ',(/ 'lev' /), 'I', 'mmr ','DPIE_O2' , gridname='physgrid')
84 0 : call addfld ('DPIE_O ',(/ 'lev' /), 'I', 'mmr ','DPIE_O' , gridname='physgrid')
85 0 : call addfld ('DPIE_N2 ',(/ 'lev' /), 'I', 'mmr ','DPIE_N2' , gridname='physgrid')
86 0 : call addfld ('DPIE_TE ',(/ 'lev' /), 'I', 'deg K ','DPIE_TE' , gridname='physgrid')
87 0 : call addfld ('DPIE_TI ',(/ 'lev' /), 'I', 'deg K ','DPIE_TI' , gridname='physgrid')
88 :
89 0 : call addfld ('PED_phys',(/ 'lev' /), 'I', 'S/m','Pedersen Conductivity' , gridname='physgrid')
90 0 : call addfld ('HAL_phys',(/ 'lev' /), 'I', 'S/m','Hall Conductivity' , gridname='physgrid')
91 :
92 0 : call addfld ('DPIE_OPMMR' ,(/ 'lev' /), 'I', 'mmr' ,'DPIE_OPMMR' , gridname='physgrid')
93 0 : call addfld ('DPIE_O2P',(/ 'lev' /), 'I', 'm^-3','DPIE_O2P(dpie input)', gridname='physgrid')
94 0 : call addfld ('DPIE_NOP',(/ 'lev' /), 'I', 'm^-3','DPIE_NOP(dpie input)', gridname='physgrid')
95 0 : call addfld ('DPIE_N2P',(/ 'lev' /), 'I', 'm^-3','DPIE_N2P(dpie input)', gridname='physgrid')
96 :
97 0 : call addfld ('WACCM_UI' ,(/ 'lev' /), 'I', 'm/s' ,'WACCM_UI (dpie output)', gridname='physgrid')
98 0 : call addfld ('WACCM_VI' ,(/ 'lev' /), 'I', 'm/s' ,'WACCM_VI (dpie output)', gridname='physgrid')
99 0 : call addfld ('WACCM_WI' ,(/ 'lev' /), 'I', 'm/s' ,'WACCM_WI (dpie output)', gridname='physgrid')
100 0 : call addfld ('WACCM_OP' ,(/ 'lev' /), 'I', 'kg/kg' ,'WACCM_OP (dpie output)', gridname='physgrid')
101 :
102 0 : call addfld ('EDYN_ADOTV1 ', (/ 'lev' /), 'I', ' ','EDYN_ADOTV1' , gridname='geo_grid')
103 0 : call addfld ('EDYN_ADOTV2 ', (/ 'lev' /), 'I', ' ','EDYN_ADOTV2' , gridname='geo_grid')
104 0 : call addfld ('EDYN_ADOTA1 ', horiz_only , 'I', ' ','EDYN_ADOTA1' , gridname='geo_grid')
105 0 : call addfld ('EDYN_ADOTA2 ', horiz_only , 'I', ' ','EDYN_ADOTA2' , gridname='geo_grid')
106 0 : call addfld ('EDYN_A1DTA2 ', horiz_only , 'I', ' ','EDYN_A1DTA2' , gridname='geo_grid')
107 :
108 0 : call addfld ('EDYN_SINI ', horiz_only , 'I', ' ','EDYN_SINI' , gridname='geo_grid')
109 0 : call addfld ('EDYN_BE3 ', horiz_only , 'I', ' ','EDYN_BE3' , gridname='geo_grid')
110 :
111 0 : call addfld ('ADOTA1_MAG', horiz_only , 'I', ' ','ADOTA1 in geo-mag coords' , gridname='gmag_grid')
112 0 : call addfld ('SINI_MAG', horiz_only , 'I', ' ','sini in geo-mag coords' , gridname='gmag_grid')
113 :
114 0 : call addfld ('OPLUS', (/ 'lev' /), 'I', 'cm^3','O+ (oplus_xport output)', gridname='geo_grid')
115 0 : call addfld ('OPtm1i',(/ 'lev' /), 'I', 'cm^3','O+ (oplus_xport output)', gridname='geo_grid')
116 0 : call addfld ('OPtm1o',(/ 'lev' /), 'I', 'cm^3','O+ (oplus_xport output)', gridname='geo_grid')
117 : endif
118 :
119 768 : end subroutine d_pie_init
120 :
121 : !-----------------------------------------------------------------------
122 8064 : subroutine d_pie_epotent( highlat_potential_model, crit_out, cols, cole, efx_phys, kev_phys, amie_in, ltr_in )
123 : use edyn_solve, only: pfrac ! NH fraction of potential (nmlonp1,nmlat0)
124 : use time_manager, only: get_curr_date
125 : use heelis, only: heelis_model
126 : use wei05sc, only: weimer05 ! driver for weimer high-lat convection model
127 : use edyn_esmf, only: edyn_esmf_update
128 : use solar_parms_data, only: solar_parms_advance
129 : use solar_wind_data, only: solar_wind_advance
130 : use solar_wind_data, only: bzimf=>solar_wind_bzimf
131 : use solar_wind_data, only: byimf=>solar_wind_byimf
132 : use solar_wind_data, only: swvel=>solar_wind_swvel
133 : use solar_wind_data, only: swden=>solar_wind_swden
134 : use edyn_mpi, only: mlat0, mlat1, mlon0, mlon1, omlon1, ntask, mytid
135 : use edyn_maggrid, only: nmlonp1, nmlat
136 : use regridder, only: regrid_mag2phys_2d
137 :
138 : ! Args:
139 : !
140 : character(len=*), intent(in) :: highlat_potential_model
141 : real(r8), intent(out) :: crit_out(2) ! critical colatitudes (degrees)
142 : integer, optional, intent(in) :: cols, cole
143 : logical, optional,intent(in) :: amie_in
144 : logical, optional,intent(in) :: ltr_in
145 :
146 : ! Prescribed energy flux
147 : real(r8), optional, intent(out) :: efx_phys(:)
148 : ! Prescribed characteristic mean energy
149 : real(r8), optional, intent(out) :: kev_phys(:)
150 :
151 : !
152 : ! local vars
153 : !
154 : logical :: amie_inputs, ltr_inputs
155 :
156 : real(r8) :: secs ! time of day in seconds
157 : integer :: iyear,imo,iday,tod ! tod is time-of-day in seconds
158 : real(r8) :: sunlon
159 :
160 : integer :: iprint
161 : integer :: j, iamie, iltr, ierr
162 : !
163 : ! AMIE fields (extra dimension added for longitude switch)
164 : !
165 16128 : real(r8) :: prescr_efxm(nmlonp1,nmlat), prescr_kevm(nmlonp1,nmlat)
166 16128 : real(r8) :: prescr_phihm(nmlonp1,nmlat)
167 :
168 8064 : call edyn_esmf_update()
169 :
170 8064 : call get_curr_date(iyear, imo,iday, tod)
171 : ! tod is integer time-of-day in seconds
172 8064 : secs = real(tod, r8)
173 :
174 : ! update solar wind data (IMF, etc.)
175 8064 : call solar_wind_advance()
176 :
177 : ! update kp -- phys timestep init happens later ...
178 8064 : call solar_parms_advance()
179 8064 : if ( mytid<ntask ) then
180 : !
181 : ! Get sun's longitude at latitudes (geographic):
182 : !
183 8064 : call sunloc(iday, secs, sunlon) ! sunlon is returned
184 : !
185 : ! Get high-latitude convection from empirical model (heelis or weimer).
186 : ! High-latitude potential phihm (edyn_solve) is defined for edynamo.
187 : !
188 8064 : if (trim(highlat_potential_model) == 'heelis') then
189 8064 : call heelis_model(sunlon) ! heelis.F90
190 0 : elseif (trim(highlat_potential_model) == 'weimer') then
191 : !
192 0 : call weimer05(byimf, bzimf, swvel, swden, sunlon)
193 0 : if (debug .and. masterproc) then
194 : write(iulog, "(a,2f8.2,a,2f8.2)") &
195 0 : 'dpie_coupling call weimer05: byimf,bzimf=', &
196 0 : byimf, bzimf, ' swvel,swden=', swvel, swden
197 : end if
198 : else
199 0 : call endrun('d_pie_epotent: Unknown highlat_potential_model')
200 : end if
201 : endif
202 :
203 8064 : amie_inputs=.false.
204 8064 : ltr_inputs=.false.
205 8064 : if (present(amie_in)) amie_inputs=amie_in
206 8064 : if (present(ltr_in)) ltr_inputs= ltr_in
207 :
208 8064 : prescribed_inputs: if (amie_inputs .or. ltr_inputs) then
209 :
210 0 : if (.not. (present(kev_phys).and.present(efx_phys)) ) then
211 0 : call endrun('d_pie_epotent: kev_phys and efx_phys must be present')
212 : end if
213 :
214 0 : iprint = 1
215 0 : if (amie_inputs) then
216 0 : if (masterproc) then
217 0 : write(iulog,*) 'Calling getamie >>> '
218 : end if
219 :
220 : call getamie(iyear, imo, iday, tod, sunlon, iprint, iamie, &
221 0 : prescr_phihm, prescr_efxm, prescr_kevm, crad)
222 :
223 0 : if (masterproc) then
224 0 : write(iulog,"('After Calling getamie >>> iamie = ', i2)") iamie
225 : end if
226 0 : prescribed_period = iamie == 1
227 : else
228 0 : if (masterproc) then
229 0 : write(iulog,*) 'Calling getltr >>> '
230 : end if
231 :
232 : call getltr(iyear, imo, iday, tod,sunlon, iprint, iltr, &
233 0 : prescr_phihm, prescr_efxm, prescr_kevm )
234 :
235 0 : if (masterproc) then
236 0 : write(iulog,"('After Calling getltr >>> iltr = ', i2)") iltr
237 : end if
238 0 : prescribed_period = iltr == 1
239 : end if
240 :
241 0 : do j = mlat0, mlat1
242 0 : call outfld('prescr_phihm',prescr_phihm(mlon0:omlon1,j),omlon1-mlon0+1,j)
243 0 : call outfld('prescr_efxm', prescr_efxm(mlon0:omlon1,j), omlon1-mlon0+1,j)
244 0 : call outfld('prescr_kevm', prescr_kevm(mlon0:omlon1,j), omlon1-mlon0+1,j)
245 : end do
246 :
247 0 : if (prescribed_period) then
248 0 : phihm = prescr_phihm
249 : end if
250 :
251 0 : call mpi_bcast(prescribed_period, 1, mpi_logical, masterprocid, mpicom, ierr)
252 :
253 0 : call regrid_mag2phys_2d(prescr_kevm(mlon0:mlon1,mlat0:mlat1), kev_phys, cols, cole)
254 0 : call regrid_mag2phys_2d(prescr_efxm(mlon0:mlon1,mlat0:mlat1), efx_phys, cols, cole)
255 :
256 0 : call outfld_phys1d( 'prescr_efxp', efx_phys )
257 0 : call outfld_phys1d( 'prescr_kevp', kev_phys )
258 :
259 : end if prescribed_inputs
260 :
261 8064 : if ( mytid<ntask ) then
262 :
263 8064 : call calc_pfrac(sunlon, pfrac) ! returns pfrac for dynamo (edyn_solve)
264 :
265 24192 : crit_out(:) = crit(:)*rtd ! degrees
266 : endif
267 8064 : end subroutine d_pie_epotent
268 :
269 : !-----------------------------------------------------------------------
270 14592 : subroutine d_pie_coupling(omega, pmid, zgi, zht, u, v, tn, &
271 7296 : sigma_ped, sigma_hall, te, ti, mbar, n2mmr, o2mmr, o1mmr, o2pmmr, &
272 7296 : nopmmr, n2pmmr, opmmr, opmmrtm1, ui, vi, wi, &
273 : rmassO2p, rmassNOp, rmassN2p, rmassOp, cols, cole, plev )
274 : !
275 : ! Call dynamo to calculate electric potential and electric field
276 : ! Note: dynamo calculates ion drifts.
277 : ! Then call oplus_xport to transport O+, which is passed back to physics.
278 : !
279 8064 : use edyn_geogrid, only: nlev
280 : use shr_const_mod, only: grav => shr_const_g ! gravitational accel. (m/s^2)
281 : use shr_const_mod, only: kboltz => shr_const_boltz ! Boltzmann const. (J/K/molecule)
282 : use time_manager, only: get_nstep, get_curr_date
283 : use edynamo, only: dynamo
284 : use edyn_mpi, only: mp_geo_halos, mp_pole_halos
285 : ! Mag grid distribution info
286 : use edyn_mpi, only: mlon0, mlon1, mlat0, mlat1, mlev0, mlev1
287 : use edyn_mpi, only: lon0, lon1, lat0, lat1, lev0, lev1, ntask, mytid
288 : use oplus, only: oplus_xport
289 : use ref_pres, only: pref_mid
290 : use regridder, only: regrid_phys2geo_3d, regrid_phys2mag_3d, regrid_geo2phys_3d
291 : use regridder, only: regrid_geo2mag_3d, regrid_geo2mag_2d
292 : use adotv_mod, only: calc_adotv
293 :
294 : !
295 : ! Args:
296 : !
297 : integer,intent(in) :: &
298 : cols, & ! First column (usually 1)
299 : cole, & ! Last column
300 : plev ! Number of levels in physics midpoint fields
301 :
302 : ! Input arguments on physics grid
303 : real(r8), intent(in) :: omega(plev, cols:cole) ! pressure velocity on midpoints (Pa/s) (i,k,j)
304 : real(r8), intent(in) :: pmid(plev, cols:cole) ! pressure at midpoints (Pa)
305 : real(r8), intent(in) :: zgi(plev, cols:cole) ! geopotential height (on interfaces) (m)
306 : real(r8), intent(in) :: zht(plev, cols:cole) ! geometric height (m) (Simple method - interfaces)
307 : real(r8), intent(in) :: u(plev, cols:cole) ! U-wind (m/s)
308 : real(r8), intent(in) :: v(plev, cols:cole) ! V-wind (m/s)
309 : real(r8), intent(in) :: tn(plev, cols:cole) ! neutral temperature (K)
310 : real(r8), intent(in) :: sigma_ped(plev, cols:cole) ! Pedersen conductivity
311 : real(r8), intent(in) :: sigma_hall(plev, cols:cole) ! Hall conductivity
312 : real(r8), intent(in) :: te(plev, cols:cole) ! electron temperature
313 : real(r8), intent(in) :: ti(plev, cols:cole) ! ion temperature
314 : real(r8), intent(in) :: mbar(plev, cols:cole) ! mean molecular weight kg/kmole
315 : real(r8), intent(in) :: n2mmr(plev, cols:cole) ! N2 mass mixing ratio (for oplus)
316 : real(r8), intent(in) :: o2mmr(plev, cols:cole) ! O2 mass mixing ratio (for oplus)
317 : real(r8), intent(in) :: o1mmr(plev, cols:cole) ! O mass mixing ratio (for oplus)
318 : real(r8), intent(in) :: o2pmmr(plev, cols:cole) ! O2+ mass mixing ratio (for oplus)
319 : real(r8), intent(in) :: nopmmr(plev, cols:cole) ! NO+ mass mixing ratio (for oplus)
320 : real(r8), intent(in) :: n2pmmr(plev, cols:cole) ! N2+ mass mixing ratio (for oplus)
321 : real(r8), intent(inout) :: opmmr(plev, cols:cole) ! O+ mass mixing ratio (oplus_xport output)
322 : real(r8), intent(inout) :: opmmrtm1(plev, cols:cole) ! O+ previous time step (oplus_xport output)
323 : real(r8), intent(inout) :: ui(plev, cols:cole) ! zonal ion drift (edynamo or empirical)
324 : real(r8), intent(inout) :: vi(plev, cols:cole) ! meridional ion drift (edynamo or empirical)
325 : real(r8), intent(inout) :: wi(plev, cols:cole) ! vertical ion drift (edynamo or empirical)
326 : real(r8), intent(in) :: rmassO2p ! O2+ molecular weight kg/kmol
327 : real(r8), intent(in) :: rmassNOp ! NO+ molecular weight kg/kmol
328 : real(r8), intent(in) :: rmassN2p ! N2+ molecular weight kg/kmol
329 : real(r8), intent(in) :: rmassOp ! O+ molecular weight kg/kmol
330 : !
331 : ! Local:
332 : !
333 : integer :: i, j, k, kk
334 : integer :: kx ! Vertical index at peak of F2 layer electron density
335 : integer :: nstep
336 : integer :: nfields ! Number of fields for multi-field calls
337 : integer :: iyear,imo,iday,tod ! tod is time-of-day in seconds
338 : integer :: isplit ! loop index
339 :
340 : real(r8) :: secs ! time of day in seconds
341 :
342 : real(r8), parameter :: small = 1.e-25_r8 ! for fields not currently available
343 :
344 14592 : real(r8) :: wn(plev, cols:cole) ! vertical velocity (from omega)
345 14592 : real(r8) :: pmid_inv(nlev) ! inverted reference pressure at midpoints (Pa)
346 :
347 : real(r8),dimension(plev, cols:cole) :: & ! ion number densities (m^3)
348 14592 : o2p, nop, n2p, op, ne, optm1
349 :
350 14592 : real(r8) :: op_in(nlev,lon0:lon1,lat0:lat1)
351 14592 : real(r8) :: optm1_in(nlev,lon0:lon1,lat0:lat1)
352 14592 : real(r8) :: zpot_in(nlev,lon0:lon1,lat0:lat1) ! geopotential height (cm)
353 14592 : real(r8) :: ui_in(nlev,lon0:lon1,lat0:lat1)
354 14592 : real(r8) :: vi_in(nlev,lon0:lon1,lat0:lat1)
355 14592 : real(r8) :: wi_in(nlev,lon0:lon1,lat0:lat1)
356 14592 : real(r8) :: wn_in(nlev,lon0:lon1,lat0:lat1) ! vertical velocity (cm/s)
357 14592 : real(r8) :: op_out(nlev,lon0:lon1,lat0:lat1)
358 14592 : real(r8) :: optm1_out(nlev,lon0:lon1,lat0:lat1)
359 :
360 14592 : real(r8),target :: halo_tn(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! neutral temperature (deg K)
361 14592 : real(r8),target :: halo_te(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! electron temperature (deg K)
362 14592 : real(r8),target :: halo_ti(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! ion temperature (deg K)
363 14592 : real(r8),target :: halo_un(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! neutral zonal wind (cm/s)
364 14592 : real(r8),target :: halo_vn(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! neutral meridional wind (cm/s)
365 14592 : real(r8),target :: halo_om(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! omega (1/s)
366 14592 : real(r8),target :: halo_o2(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! o2 (mmr)
367 14592 : real(r8),target :: halo_o1(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! o (mmr)
368 14592 : real(r8),target :: halo_n2(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! n2 (mmr)
369 14592 : real(r8),target :: halo_mbar(nlev,lon0-2:lon1+2,lat0-2:lat1+2) ! mean molecular weight
370 7296 : real(r8), allocatable :: polesign(:)
371 : !
372 14592 : real(r8) :: nmf2(cols:cole) ! Electron number density at F2 peak (m-3 converted to cm-3)
373 14592 : real(r8) :: hmf2(cols:cole) ! Height of electron number density F2 peak (m converted to km)
374 : real(r8) :: &
375 : height(3), & ! Surrounding heights when locating electron density F2 peak
376 : nde(3) ! Surround densities when locating electron density F2 peak
377 : real(r8) h12,h22,h32,deltx,atx,ax,btx,bx,ctx,cx ! Variables used for weighting when locating F2 peak
378 : !
379 : logical :: do_integrals
380 : !
381 : ! Pointers for multiple-field calls:
382 7296 : type(array_ptr_type),allocatable :: ptrs(:)
383 :
384 : character(len=*), parameter :: subname = 'd_pie_coupling'
385 :
386 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) :: & ! 3d fields on mag grid
387 14592 : ped_mag, & ! pedersen conductivity on magnetic grid
388 14592 : hal_mag, & ! hall conductivity on magnetic grid
389 14592 : zpot_mag ! geopotential on magnetic grid
390 :
391 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) :: & ! 3d fields on mag grid
392 14592 : ped_mag_in, & ! pedersen conductivity on magnetic grid
393 14592 : hal_mag_in, & ! hall conductivity on magnetic grid
394 14592 : zpot_mag_in ! geopotential on magnetic grid
395 :
396 : real(r8), dimension(lon0:lon1,lat0:lat1,lev0:lev1) :: & ! 3d fields on geo grid
397 14592 : zpot_geo, & ! geopotential on magnetic grid
398 14592 : tn_geo, &
399 14592 : te_geo, &
400 14592 : ti_geo, &
401 14592 : un_geo, &
402 14592 : vn_geo, &
403 14592 : wn_geo, &
404 14592 : ui_geo, &
405 14592 : vi_geo, &
406 14592 : wi_geo, &
407 14592 : omega_geo, &
408 14592 : o2_geo, &
409 14592 : o_geo, &
410 14592 : n2_geo, &
411 14592 : op_geo, &
412 14592 : optm1_geo, &
413 14592 : pmid_geo, &
414 14592 : mbar_geo
415 :
416 : real(r8), dimension(lon0:lon1,lat0:lat1,lev0:lev1) :: &
417 14592 : adotv1_in, adotv2_in
418 : real(r8), dimension(lon0:lon1,lat0:lat1) :: &
419 21888 : adota1_in, adota2_in, a1dta2_in, be3_in, sini_in
420 :
421 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1) :: &
422 14592 : adotv1_mag, adotv2_mag
423 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1) :: &
424 21888 : adota1_mag, adota2_mag, a1dta2_mag, be3_mag, sini_mag
425 :
426 7296 : call t_startf(subname)
427 :
428 7296 : if (debug .and. masterproc) then
429 :
430 0 : nstep = get_nstep()
431 0 : call get_curr_date(iyear, imo, iday, tod)
432 0 : secs = real(tod, r8)
433 :
434 : write(iulog,"(3a,i8,a,3i5,a,f6.2)") &
435 0 : 'Enter ',subname,': nstep = ',nstep, ', iyear,imo,iday = ', &
436 0 : iyear, imo, iday, ' ut (hrs) = ', secs/3600._r8
437 :
438 0 : write(iulog,"(a,': nspltop = ',i3)") subname, nspltop
439 : end if
440 :
441 : !---------------------------------------------------------------
442 : ! Calculate vertical neutral wind velocity wn(i,j,k).
443 : ! (omega (Pa/s), tn (K), and mbar (kg/kmole) are inputs, grav is m/s^2)
444 : !---------------------------------------------------------------
445 7296 : call calc_wn(tn, omega, pmid, mbar, grav, wn, cols, cole, nlev)
446 :
447 : !---------------------------------------------------------------
448 : ! Convert from mmr to number densities (m^3):
449 : !---------------------------------------------------------------
450 269952 : do i = cols, cole
451 34415232 : do k = 1, nlev
452 : ! O2+, NO+, N2+, O+:
453 68290560 : o2p(k, i) = o2pmmr(k, i) * mbar(k, i) / rmassO2p * &
454 68290560 : pmid(k,i) / (kboltz * tn(k, i))
455 : nop(k, i) = nopmmr(k, i) * mbar(k, i) / rmassNOp * &
456 34145280 : pmid(k,i) / (kboltz * tn(k, i))
457 : n2p(k, i) = n2pmmr(k, i) * mbar(k, i) / rmassN2p * &
458 34145280 : pmid(k,i) / (kboltz * tn(k, i))
459 : op(k, i) = opmmr(k, i) * mbar(k, i) / rmassOp * &
460 34145280 : pmid(k,i) / (kboltz * tn(k, i))
461 : optm1(k, i) = opmmrtm1(k, i) * mbar(k, i) / rmassOp * &
462 34145280 : pmid(k,i) / (kboltz * tn(k, i))
463 34407936 : ne(k, i) = o2p(k,i)+nop(k,i)+n2p(k,i)+op(k,i)
464 : end do
465 : end do ! k=1,nlev
466 :
467 7296 : if (debug_hist) then
468 0 : call outfld_phys('DPIE_TN',tn)
469 0 : call outfld_phys('DPIE_UN',u* 100._r8)
470 0 : call outfld_phys('DPIE_VN',v* 100._r8)
471 0 : call outfld_phys('DPIE_WN',wn* 100._r8)
472 0 : call outfld_phys('DPIE_ZHT',zht* 100._r8)
473 0 : call outfld_phys('DPIE_ZGI',zgi* 100._r8)
474 0 : call outfld_phys('DPIE_MBAR',mbar)
475 :
476 0 : call outfld_phys('DPIE_N2',n2mmr)
477 0 : call outfld_phys('DPIE_O2',o2mmr)
478 0 : call outfld_phys('DPIE_O',o1mmr)
479 :
480 0 : call outfld_phys('DPIE_OMEGA',omega)
481 0 : call outfld_phys('DPIE_OM',-omega/pmid)
482 :
483 0 : call outfld_phys('DPIE_TE',te)
484 0 : call outfld_phys('DPIE_TI',ti)
485 :
486 0 : call outfld_phys('DPIE_O2P',o2p)
487 0 : call outfld_phys('DPIE_NOP',nop)
488 0 : call outfld_phys('DPIE_N2P',n2p)
489 : endif
490 34415232 : call outfld_phys('EDens',ne/1.E6_r8)
491 34415232 : call outfld_phys('OpDens',op/1.E6_r8)
492 :
493 : !-------------------------------------------------------------------------
494 : ! Derive diagnostics nmF2 and hmF2 for output based on TIE-GCM algorithm
495 : !-------------------------------------------------------------------------
496 7296 : if (hist_fld_active('HMF2') .or. hist_fld_active('NMF2')) then
497 0 : iloop: do i = cols, cole
498 0 : kx = 0
499 0 : kloop: do k= 2, nlev
500 0 : if (ne(k,i) >= ne(k-1,i) .and. ne(k,i) >= ne(k+1,i)) then
501 : kx = k
502 : exit kloop
503 : end if
504 : end do kloop
505 :
506 0 : if (kx==0) then
507 0 : hmf2(i) = fillvalue
508 0 : nmf2(i) = fillvalue
509 0 : exit iloop
510 : end if
511 :
512 0 : height = (/zht(kx+1,i),zht(kx,i),zht(kx-1,i)/)
513 0 : nde = (/ne(kx+1,i),ne(kx,i),ne(kx-1,i)/)
514 :
515 0 : h12 = height(1)*height(1)
516 0 : h22 = height(2)*height(2)
517 0 : h32 = height(3)*height(3)
518 :
519 0 : deltx=h12*height(2)+h22*height(3)+h32*height(1)-h32*height(2)-h12*height(3)-h22*height(1)
520 :
521 0 : atx=nde(1)*height(2)+nde(2)*height(3)+nde(3)*height(1)-height(2)*nde(3)-height(3)*nde(1)-height(1)*nde(2)
522 0 : ax=atx/deltx
523 :
524 0 : btx=h12*nde(2)+h22*nde(3)+h32*nde(1)-h32*nde(2)-h12*nde(3)-h22*nde(1)
525 0 : bx=btx/deltx
526 : ctx=h12*height(2)*nde(3)+h22*height(3)*nde(1)+h32*height(1)*nde(2)-h32*height(2)*nde(1)- &
527 0 : h12*height(3)*nde(2)-h22*height(1)*nde(3)
528 0 : cx=ctx/deltx
529 :
530 0 : hmf2(i)=-(bx/(2._r8*ax)) * 1.E-03_r8
531 0 : nmf2(i)=-((bx*bx-4._r8*ax*cx)/(4._r8*ax)) * 1.E-06_r8
532 :
533 : end do iloop ! i=cols, cole
534 :
535 0 : call outfld_phys1d('HMF2',hmf2)
536 0 : call outfld_phys1d('NMF2',nmf2)
537 : end if
538 7296 : if (debug_hist) then
539 0 : call outfld_phys('DPIE_OPMMR', opmmr)
540 0 : call outfld_phys('PED_phys', sigma_ped )
541 0 : call outfld_phys('HAL_phys', sigma_hall )
542 : endif
543 7296 : if (ionos_edyn_active .or. ionos_oplus_xport) then
544 :
545 7296 : call regrid_phys2geo_3d( zgi,zpot_geo, plev, cols, cole )
546 7296 : call regrid_phys2geo_3d( u, un_geo, plev, cols, cole )
547 7296 : call regrid_phys2geo_3d( v, vn_geo, plev, cols, cole )
548 7296 : call regrid_phys2geo_3d( wn,wn_geo, plev, cols, cole )
549 7296 : call regrid_phys2geo_3d( ui, ui_geo, plev, cols, cole )
550 7296 : call regrid_phys2geo_3d( vi, vi_geo, plev, cols, cole )
551 7296 : call regrid_phys2geo_3d( wi, wi_geo, plev, cols, cole )
552 :
553 955776 : do k = 1, nlev
554 948480 : kk = nlev-k+1
555 3801216 : do j = lat0, lat1
556 37939200 : do i = lon0, lon1
557 34145280 : zpot_in(kk,i,j) = zpot_geo(i,j,k) * 100._r8 ! m -> cm
558 34145280 : halo_un(kk,i,j) = un_geo(i,j,k) * 100._r8 ! m/s -> cm/s
559 34145280 : halo_vn(kk,i,j) = vn_geo(i,j,k) * 100._r8 ! m/s -> cm/s
560 34145280 : wn_in(kk,i,j) = wn_geo(i,j,k) * 100._r8 ! m/s -> cm/s
561 34145280 : ui_in(kk,i,j) = ui_geo(i,j,k) * 100._r8 ! zonal ion drift (m/s -> cm/s)
562 34145280 : vi_in(kk,i,j) = vi_geo(i,j,k) * 100._r8 ! meridional ion drift (m/s -> cm/s)
563 36990720 : wi_in(kk,i,j) = wi_geo(i,j,k) * 100._r8 ! vertical ion drift (m/s -> cm/s)
564 : end do
565 : end do
566 : end do
567 :
568 : end if
569 :
570 : !
571 : !
572 : ! Call electrodynamo (edynamo.F90)
573 : ! If using time3d conductances, tell dynamo to *not* do fieldline
574 : ! integrations (i.e., do_integrals == false). In this case, edynamo
575 : ! conductances zigmxx,rim1,2 from time3d will be set by subroutine
576 : ! transform_glbin in time3d module.
577 : !
578 7296 : do_integrals = .true.
579 : !
580 : ! If ionos_edyn_active=false, then empirical ion drifts were passed in from physics,
581 : ! otherwise dynamo calculates them here, and they will be passed to physics.
582 : !
583 7296 : if (ionos_edyn_active) then
584 :
585 7296 : call t_startf('dpie_ionos_dynamo')
586 :
587 0 : call calc_adotv( zpot_in(lev0:lev1,lon0:lon1,lat0:lat1), &
588 0 : halo_un(lev0:lev1,lon0:lon1,lat0:lat1), &
589 : halo_vn(lev0:lev1,lon0:lon1,lat0:lat1), &
590 : wn_in(lev0:lev1,lon0:lon1,lat0:lat1), &
591 : adotv1_in, adotv2_in, adota1_in, adota2_in, &
592 68866944 : a1dta2_in, be3_in, sini_in, lev0, lev1, lon0, lon1, lat0, lat1)
593 :
594 7296 : call regrid_geo2mag_3d( adotv1_in, adotv1_mag )
595 7296 : call regrid_geo2mag_3d( adotv2_in, adotv2_mag )
596 7296 : if (debug_hist) then
597 0 : call outfld_geo('EDYN_ADOTV1', adotv1_in(:,:,lev1:lev0:-1) )
598 0 : call outfld_geo('EDYN_ADOTV2', adotv2_in(:,:,lev1:lev0:-1) )
599 :
600 0 : call outfld_geo2d( 'EDYN_ADOTA1', adota1_in )
601 0 : call outfld_geo2d( 'EDYN_ADOTA2', adota2_in )
602 0 : call outfld_geo2d( 'EDYN_A1DTA2', a1dta2_in )
603 0 : call outfld_geo2d( 'EDYN_BE3' , be3_in )
604 0 : call outfld_geo2d( 'EDYN_SINI', sini_in )
605 : endif
606 7296 : call regrid_geo2mag_2d( adota1_in, adota1_mag )
607 7296 : call regrid_geo2mag_2d( adota2_in, adota2_mag )
608 7296 : call regrid_geo2mag_2d( a1dta2_in, a1dta2_mag )
609 7296 : call regrid_geo2mag_2d( be3_in, be3_mag )
610 7296 : call regrid_geo2mag_2d( sini_in, sini_mag )
611 7296 : if (debug_hist) then
612 0 : call outfld_mag2d('ADOTA1_MAG', adota1_mag )
613 0 : call outfld_mag2d('SINI_MAG', sini_mag )
614 : endif
615 7296 : call regrid_phys2mag_3d( sigma_ped, ped_mag, plev, cols, cole )
616 7296 : call regrid_phys2mag_3d( sigma_hall, hal_mag, plev, cols, cole )
617 7296 : call regrid_phys2mag_3d( zgi, zpot_mag, plev, cols, cole )
618 :
619 7296 : if (mytid<ntask) then
620 23237646 : zpot_mag_in(:,:,mlev0:mlev1) = zpot_mag(:,:,mlev1:mlev0:-1) * 100._r8 ! m -> cm
621 23237646 : ped_mag_in(:,:,mlev0:mlev1) = ped_mag(:,:,mlev1:mlev0:-1)
622 23237646 : hal_mag_in(:,:,mlev0:mlev1) = hal_mag(:,:,mlev1:mlev0:-1)
623 :
624 : call dynamo( zpot_mag_in, ped_mag_in, hal_mag_in, adotv1_mag, adotv2_mag, adota1_mag, &
625 : adota2_mag, a1dta2_mag, be3_mag, sini_mag, &
626 : zpot_in, ui_in, vi_in, wi_in, &
627 7296 : lon0,lon1, lat0,lat1, lev0,lev1, do_integrals )
628 : endif
629 :
630 7296 : call t_stopf ('dpie_ionos_dynamo')
631 :
632 : else
633 0 : if (debug .and. masterproc) then
634 0 : write(iulog,"('dpie_coupling (dynamo NOT called): nstep=',i8)") nstep
635 : write(iulog,"(' empirical ExB ui min,max (cm/s)=',2es12.4)") &
636 0 : minval(ui),maxval(ui)
637 : write(iulog,"(' empirical ExB vi min,max (cm/s)=',2es12.4)") &
638 0 : minval(vi),maxval(vi)
639 : write(iulog,"(' empirical ExB wi min,max (cm/s)=',2es12.4)") &
640 0 : minval(wi),maxval(wi)
641 : end if
642 : end if
643 :
644 : !
645 : ! Call O+ transport routine. Now all inputs to oplus_xport should be in
646 : ! tiegcm-format wrt longitude (-180->180), vertical (bot2top), and units (CGS).
647 : ! (Composition is mmr, ne is cm^3, winds are cm/s)
648 : ! Output op_out and opnm_out will be in cm^3, converted to mmr below.
649 : !
650 7296 : if (ionos_oplus_xport) then
651 955776 : pmid_inv(1:nlev) = pref_mid(nlev:1:-1) ! invert ref pressure (Pa) as in tiegcm
652 :
653 :
654 : !
655 : ! Transport O+ (all args in 'TIEGCM format')
656 : ! Subcycle oplus_xport nspltop times.
657 : !
658 7296 : if (debug .and. masterproc) then
659 : write(iulog,"(a,i8,a,i3)") &
660 0 : 'dpie_coupling before subcycling oplus_xport: nstep = ', &
661 0 : nstep, ' nspltop = ', nspltop
662 : end if
663 :
664 7296 : call regrid_phys2geo_3d( tn, tn_geo, plev, cols, cole )
665 7296 : call regrid_phys2geo_3d( te, te_geo, plev, cols, cole )
666 7296 : call regrid_phys2geo_3d( ti, ti_geo, plev, cols, cole )
667 7296 : call regrid_phys2geo_3d( omega, omega_geo, plev, cols, cole )
668 7296 : call regrid_phys2geo_3d( o2mmr, o2_geo, plev, cols, cole )
669 7296 : call regrid_phys2geo_3d( n2mmr, n2_geo, plev, cols, cole )
670 7296 : call regrid_phys2geo_3d( o1mmr, o_geo, plev, cols, cole )
671 7296 : call regrid_phys2geo_3d( op, op_geo, plev, cols, cole )
672 7296 : call regrid_phys2geo_3d( optm1, optm1_geo, plev, cols, cole )
673 7296 : call regrid_phys2geo_3d( pmid, pmid_geo, plev, cols, cole )
674 7296 : call regrid_phys2geo_3d( mbar, mbar_geo, plev, cols, cole )
675 :
676 7296 : if (mytid<ntask) then
677 :
678 7296 : call t_startf('dpie_halo')
679 955776 : do k = 1, nlev
680 948480 : kk = nlev-k+1
681 3801216 : do j = lat0, lat1
682 37939200 : do i = lon0, lon1
683 34145280 : halo_tn(kk,i,j) = tn_geo(i,j,k)
684 34145280 : halo_te(kk,i,j) = te_geo(i,j,k)
685 34145280 : halo_ti(kk,i,j) = ti_geo(i,j,k)
686 34145280 : halo_om(kk,i,j) = -omega_geo(i,j,k) / pmid_geo(i,j,k) ! Pa/s -> 1/s
687 34145280 : halo_o2(kk,i,j) = o2_geo(i,j,k)
688 34145280 : halo_o1(kk,i,j) = o_geo(i,j,k)
689 34145280 : halo_n2(kk,i,j) = n2_geo(i,j,k)
690 34145280 : halo_mbar(kk,i,j) = mbar_geo(i,j,k)
691 34145280 : op_in(kk,i,j) = op_geo(i,j,k) / 1.e6_r8 ! m^3 -> cm^3
692 36990720 : optm1_in(kk,i,j) = optm1_geo(i,j,k) / 1.e6_r8 ! m^3 -> cm^3
693 : end do
694 : end do
695 : end do
696 : !
697 : ! Define halo points on inputs:
698 : ! WACCM has global longitude values at the poles (j=1,j=nlev)
699 : ! (they are constant for most, except the winds.)
700 : !
701 : ! Set two halo points in lat,lon:
702 : !
703 7296 : nfields = 10
704 7296 : allocate(ptrs(nfields), polesign(nfields))
705 7296 : ptrs(1)%ptr => halo_tn ; ptrs(2)%ptr => halo_te ; ptrs(3)%ptr => halo_ti
706 7296 : ptrs(4)%ptr => halo_un ; ptrs(5)%ptr => halo_vn ; ptrs(6)%ptr => halo_om
707 7296 : ptrs(7)%ptr => halo_o2 ; ptrs(8)%ptr => halo_o1 ; ptrs(9)%ptr => halo_n2
708 7296 : ptrs(10)%ptr => halo_mbar
709 :
710 80256 : polesign = 1._r8
711 21888 : polesign(4:5) = -1._r8 ! un,vn
712 :
713 7296 : call mp_geo_halos(ptrs,1,nlev,lon0,lon1,lat0,lat1,nfields)
714 : !
715 : ! Set latitude halo points over the poles (this does not change the poles).
716 : ! (the 2nd halo over the poles will not actually be used (assuming lat loops
717 : ! are lat=2,plat-1), because jp1,jm1 will be the pole itself, and jp2,jm2
718 : ! will be the first halo over the pole)
719 : !
720 7296 : call mp_pole_halos(ptrs,1,nlev,lon0,lon1,lat0,lat1,nfields,polesign)
721 7296 : deallocate(ptrs,polesign)
722 7296 : call t_stopf('dpie_halo')
723 7296 : if (debug_hist) then
724 0 : call outfld_geokij( 'OPtm1i',optm1_in, lev0,lev1, lon0,lon1, lat0,lat1 )
725 : endif
726 7296 : call t_startf('dpie_oplus_xport')
727 43776 : do isplit = 1, nspltop
728 :
729 36480 : if (isplit > 1) then
730 137748480 : op_in = op_out
731 137748480 : optm1_in = optm1_out
732 : end if
733 :
734 : call oplus_xport(halo_tn, halo_te, halo_ti, halo_un, halo_vn, halo_om, &
735 : zpot_in, halo_o2, halo_o1, halo_n2, op_in, optm1_in, &
736 : halo_mbar, ui_in, vi_in, wi_in, pmid_inv, &
737 : op_out, optm1_out, &
738 43776 : lon0, lon1, lat0, lat1, nspltop, isplit)
739 :
740 : end do ! isplit=1,nspltop
741 7296 : call t_stopf ('dpie_oplus_xport')
742 7296 : if (debug.and.masterproc) then
743 : write(iulog,"('dpie_coupling after subcycling oplus_xport: nstep=',i8,' nspltop=',i3)") &
744 0 : nstep,nspltop
745 0 : write(iulog,"(' op_out min,max (cm^3)=',2es12.4)") minval(op_out) ,maxval(op_out)
746 0 : write(iulog,"(' optm1_out min,max (cm^3)=',2es12.4)") minval(optm1_out),maxval(optm1_out)
747 : end if
748 7296 : if (debug_hist) then
749 0 : call outfld_geokij( 'OPLUS', op_out, lev0,lev1, lon0,lon1, lat0,lat1 )
750 0 : call outfld_geokij( 'OPtm1o',optm1_out, lev0,lev1, lon0,lon1, lat0,lat1 )
751 : endif
752 : !
753 : ! Pass new O+ for current and previous time step back to physics (convert from cm^3 to m^3 and back to mmr).
754 : !
755 955776 : do k=1,nlev
756 948480 : kk = nlev-k+1
757 3801216 : do j = lat0,lat1
758 37939200 : do i = lon0,lon1
759 136581120 : op_geo(i,j,k) = op_out(kk,i,j)*1.e6_r8 * rmassOp / mbar_geo(i,j,k) * &
760 136581120 : (kboltz * tn_geo(i,j,k)) / pmid_geo(i,j,k)
761 : optm1_geo(i,j,k) = optm1_out(kk,i,j)*1.e6_r8 * rmassOp / mbar_geo(i,j,k) * &
762 34145280 : (kboltz * tn_geo(i,j,k)) / pmid_geo(i,j,k)
763 34145280 : ui_geo(i,j,k) = ui_in(kk,i,j)/100._r8 ! cm/s -> m/s
764 34145280 : vi_geo(i,j,k) = vi_in(kk,i,j)/100._r8 ! cm/s -> m/s
765 36990720 : wi_geo(i,j,k) = wi_in(kk,i,j)/100._r8 ! cm/s -> m/s
766 : end do
767 : end do
768 : end do
769 :
770 : endif
771 :
772 7296 : call regrid_geo2phys_3d( op_geo, opmmr, plev, cols, cole )
773 7296 : call regrid_geo2phys_3d( optm1_geo, opmmrtm1, plev, cols, cole )
774 7296 : call regrid_geo2phys_3d( ui_geo, ui, plev, cols, cole )
775 7296 : call regrid_geo2phys_3d( vi_geo, vi, plev, cols, cole )
776 7296 : call regrid_geo2phys_3d( wi_geo, wi, plev, cols, cole )
777 :
778 : end if ! ionos_oplus_xport
779 :
780 7296 : if (debug_hist) then
781 0 : call outfld_phys('WACCM_UI',ui)
782 0 : call outfld_phys('WACCM_VI',vi)
783 0 : call outfld_phys('WACCM_WI',wi)
784 0 : call outfld_phys('WACCM_OP',opmmr)
785 : endif
786 7296 : call t_stopf('d_pie_coupling')
787 :
788 7296 : end subroutine d_pie_coupling
789 : !-----------------------------------------------------------------------
790 7296 : subroutine calc_wn(tn,omega,pmid,mbar,grav,wn,cols,cole,nlev)
791 7296 : use shr_const_mod,only : shr_const_rgas ! Universal gas constant
792 : !
793 : ! Calculate neutral vertical wind on midpoints (m/s)
794 : !
795 : ! Inputs:
796 : integer,intent(in) :: cols, cole, nlev
797 : real(r8),dimension(nlev, cols:cole),intent(in) :: &
798 : tn, & ! neutral temperature (deg K)
799 : omega,& ! pressure velocity (Pa/s)
800 : mbar ! mean molecular weight
801 : real(r8),dimension(nlev,cols:cole),intent(in) :: &
802 : pmid ! pressure at midpoints (Pa)
803 : real(r8),intent(in) :: grav ! m/s^2
804 : !
805 : ! Output:
806 : real(r8),intent(out) :: wn(nlev, cols:cole) ! vertical velocity output (m/s)
807 : !
808 : ! Local:
809 : integer :: i,k
810 14592 : real(r8) :: scheight(nlev, cols:cole) ! dimensioned for vectorization
811 :
812 269952 : do i = cols, cole
813 34415232 : do k = 1, nlev
814 34145280 : scheight(k,i) = shr_const_rgas*tn(k,i)/(mbar(k,i)*grav)
815 34407936 : wn(k,i) = -omega(k,i)*scheight(k,i)/pmid(k,i)
816 : end do
817 : end do
818 7296 : end subroutine calc_wn
819 : !-----------------------------------------------------------------------
820 8064 : subroutine calc_pfrac(sunlon,pfrac)
821 : !
822 : ! Calculate pfrac fractional presence of dynamo equation using critical
823 : ! convection colatitudes crit(2).
824 : !
825 : use edyn_maggrid ,only: nmlonp1,ylonm,ylatm
826 : use edyn_solve ,only: nmlat0
827 : use aurora_params ,only: offc, dskofc, theta0, aurora_params_set
828 :
829 : implicit none
830 : !
831 : ! Args:
832 : real(r8),intent(in) :: sunlon ! Sun's longitude in dipole coordinates
833 : !
834 : ! Output: fractional presence of dynamo equation using critical colatitudes
835 : !
836 : real(r8),intent(out) :: pfrac(nmlonp1,nmlat0) ! NH fraction of potential
837 : !
838 : ! Local:
839 : integer :: j,i
840 16128 : real(r8),dimension(nmlonp1,nmlat0) :: colatc
841 : real(r8) :: sinlat,coslat,aslonc,ofdc,cosofc,sinofc,crit1deg
842 :
843 8064 : if (.not. crit_user_set) then
844 8064 : if (prescribed_period) then
845 0 : crit(:) = amie_default_crit(:)*dtr
846 : else
847 8064 : crit1deg = max(15._r8,0.5_r8*(theta0(1)+theta0(2))*rtd + 5._r8)
848 8064 : crit1deg = min(30._r8,crit1deg)
849 : !
850 : ! Critical colatitudes:
851 8064 : crit(1) = crit1deg*dtr
852 8064 : crit(2) = crit(1) + 15._r8*dtr
853 : end if
854 : end if
855 :
856 8064 : if (.not.aurora_params_set) then
857 0 : offc(:) = 1._r8*dtr
858 0 : dskofc(:) = 0._r8
859 : end if
860 :
861 : !
862 : ! offc(2), dskofc(2) are for northern hemisphere aurora
863 : !
864 8064 : ofdc = sqrt(offc(2)**2 + dskofc(2)**2)
865 8064 : cosofc = cos(ofdc)
866 8064 : sinofc = sin(ofdc)
867 8064 : aslonc = asin(dskofc(2)/ofdc)
868 : !
869 : ! Define colatc with northern convection circle coordinates
870 : !
871 403200 : do j=1,nmlat0
872 395136 : sinlat = sin(abs(ylatm(j+nmlat0-1)))
873 395136 : coslat = cos( ylatm(j+nmlat0-1))
874 32401152 : do i=1,nmlonp1
875 32006016 : colatc(i,j) = cos(ylonm(i)-sunlon+aslonc)
876 32401152 : colatc(i,j) = acos(cosofc*sinlat-sinofc*coslat*colatc(i,j))
877 : end do ! i=1,nmlonp1
878 : !
879 : ! Calculate fractional presence of dynamo equation at each northern
880 : ! hemisphere geomagnetic grid point. Output in pfrac(nmlonp1,nmlat0)
881 : !
882 32409216 : do i = 1 , nmlonp1
883 32006016 : pfrac(i,j) = (colatc(i,j)-crit(1)) / (crit(2)-crit(1))
884 32006016 : if (pfrac(i,j) < 0._r8) then
885 6633216 : pfrac(i,j) = 0._r8
886 : end if
887 32401152 : if (pfrac(i,j) >= 1._r8) then
888 21555072 : pfrac(i,j) = 1._r8
889 : end if
890 : end do ! i=1,nmlonp1
891 : end do ! j=1,nmlat0
892 : !
893 8064 : end subroutine calc_pfrac
894 : !-----------------------------------------------------------------------
895 8064 : subroutine sunloc(iday, secs, sunlon)
896 : !
897 : ! Given day of year and ut, return sun's longitude in dipole coordinates
898 : ! in sunlon
899 : !
900 : use getapex, only: alonm ! (nlonp1,0:nlatp1)
901 : use edyn_geogrid, only: nlon, nlat, dphi, dlamda
902 : use edyn_params, only: pi
903 : !
904 : ! Args:
905 : integer,intent(in) :: iday ! day of year
906 : real(r8),intent(in) :: secs ! ut in seconds
907 : real(r8),intent(out) :: sunlon ! output
908 : !
909 : ! Local:
910 : integer :: j, i, ii, isun, jsun
911 : real(r8) :: glats, glons, pisun, pjsun, sndlons, csdlons
912 16128 : real(r8) :: rlonm(nlon+4, nlat) ! (nlon+4,nlat)
913 : real(r8) :: r8_isun, r8_jsun
914 :
915 : !
916 : ! Sun's geographic coordinates:
917 8064 : glats = asin(.398749_r8*sin(2._r8 * pi * real(iday-80, r8) / 365._r8))
918 8064 : glons = pi * (1._r8 - (2._r8 * secs / 86400._r8))
919 :
920 782208 : do j = 1, nlat
921 112250880 : do i = 1, nlon
922 111476736 : ii = i + 2
923 112250880 : rlonm(ii, j) = alonm(i, j)
924 : end do
925 2330496 : do i = 1, 2
926 1548288 : rlonm(i, j) = rlonm(i+nlon, j)
927 2322432 : rlonm(i+nlon+2, j) = rlonm(i+2, j)
928 : end do
929 : end do
930 :
931 8064 : pisun = ((glons + pi) / dlamda) + 1._r8
932 8064 : pjsun = ((glats + (.5_r8 * (pi - dphi))) / dphi) + 1._r8
933 8064 : isun = int(pisun)
934 8064 : jsun = int(pjsun)
935 8064 : r8_isun = real(isun, r8)
936 8064 : r8_jsun = real(jsun, r8)
937 8064 : pisun = pisun - r8_isun
938 8064 : pjsun = pjsun - r8_jsun
939 :
940 16128 : sndlons = (1._r8-pisun) * (1._r8-pjsun) * sin(rlonm(isun+2, jsun)) + &
941 8064 : pisun*(1._r8-pjsun) * sin(rlonm(isun+3,jsun)) + &
942 8064 : pisun*pjsun * sin(rlonm(isun+3,jsun+1)) + &
943 32256 : (1._r8-pisun)*pjsun * sin(rlonm(isun+2,jsun+1))
944 : csdlons = (1._r8-pisun) * (1._r8-pjsun) * cos(rlonm(isun+2,jsun)) + &
945 : pisun*(1._r8-pjsun) * cos(rlonm(isun+3,jsun))+ &
946 : pisun*pjsun * cos(rlonm(isun+3,jsun+1))+ &
947 8064 : (1._r8-pisun)*pjsun * cos(rlonm(isun+2,jsun+1))
948 8064 : sunlon = atan2(sndlons, csdlons)
949 :
950 8064 : end subroutine sunloc
951 :
952 : !-----------------------------------------------------------------------
953 : !-----------------------------------------------------------------------
954 0 : subroutine outfld_phys1d( fldname, array )
955 : use ppgrid, only: pcols, begchunk, endchunk
956 : use phys_grid,only: get_ncols_p
957 :
958 : character(len=*), intent(in) :: fldname
959 : real(r8), intent(in) :: array(:)
960 :
961 : integer :: i,j, lchnk,ncol
962 : real(r8) :: tmparr(pcols)
963 :
964 0 : if (hist_fld_active(fldname)) then
965 0 : j = 0
966 0 : do lchnk = begchunk, endchunk
967 0 : ncol = get_ncols_p(lchnk)
968 0 : do i = 1, ncol
969 0 : j = j + 1
970 0 : tmparr(i) = array(j)
971 : enddo
972 0 : call outfld(fldname,tmparr(:ncol),ncol,lchnk)
973 : enddo
974 : end if
975 :
976 0 : end subroutine outfld_phys1d
977 : !-----------------------------------------------------------------------
978 14592 : subroutine outfld_phys( fldname, array )
979 0 : use ppgrid, only: pcols, pver, begchunk, endchunk
980 : use phys_grid,only: get_ncols_p
981 :
982 : character(len=*), intent(in) :: fldname
983 : real(r8), intent(in) :: array(:,:)
984 :
985 : integer :: i,j,k, lchnk,ncol
986 : real(r8) :: tmparr(pcols, pver)
987 :
988 14592 : if (hist_fld_active(fldname)) then
989 14592 : j = 0
990 58368 : do lchnk = begchunk, endchunk
991 43776 : ncol = get_ncols_p(lchnk)
992 569088 : do i = 1, ncol
993 525312 : j = j + 1
994 68859648 : do k = 1, pver
995 68815872 : tmparr(i,k) = array(k,j)
996 : enddo
997 : enddo
998 74039808 : call outfld(fldname,tmparr(:ncol,:),ncol,lchnk)
999 : enddo
1000 : end if
1001 :
1002 14592 : end subroutine outfld_phys
1003 : !-----------------------------------------------------------------------
1004 0 : subroutine outfld_geokij( name, array, ilev0,ilev1, ilon0,ilon1, ilat0,ilat1 )
1005 :
1006 : character(len=*), intent(in) :: name
1007 : integer, intent(in) :: ilev0,ilev1, ilon0,ilon1, ilat0,ilat1
1008 : real(r8), intent(in) :: array(ilev0:ilev1, ilon0:ilon1, ilat0:ilat1)
1009 :
1010 : integer :: j,k
1011 0 : real(r8) :: tmpout(ilon0:ilon1,ilev0:ilev1)
1012 :
1013 0 : do j = ilat0,ilat1
1014 0 : do k = ilev0,ilev1
1015 0 : tmpout(ilon0:ilon1,k) = array(ilev1-k+1,ilon0:ilon1,j)
1016 : end do
1017 0 : call outfld( name, tmpout, ilon1-ilon0+1, j )
1018 : end do
1019 14592 : end subroutine outfld_geokij
1020 : !-----------------------------------------------------------------------
1021 0 : subroutine outfld_geo( fldname, array )
1022 : use edyn_mpi, only: lon0, lon1, lat0, lat1, lev0, lev1
1023 :
1024 : character(len=*), intent(in) :: fldname
1025 : real(r8), intent(in) :: array(lon0:lon1,lat0:lat1,lev0:lev1)
1026 :
1027 : integer :: j
1028 :
1029 0 : do j = lat0,lat1
1030 0 : call outfld( fldname, array(lon0:lon1,j,lev0:lev1), lon1-lon0+1, j )
1031 : end do
1032 :
1033 0 : end subroutine outfld_geo
1034 : !-----------------------------------------------------------------------
1035 0 : subroutine outfld_geo2d( fldname, array )
1036 : use edyn_mpi, only: lon0, lon1, lat0, lat1
1037 :
1038 : character(len=*), intent(in) :: fldname
1039 : real(r8), intent(in) :: array(lon0:lon1,lat0:lat1)
1040 :
1041 : integer :: j
1042 :
1043 0 : do j = lat0,lat1
1044 0 : call outfld( fldname, array(lon0:lon1,j), lon1-lon0+1, j )
1045 : end do
1046 :
1047 0 : end subroutine outfld_geo2d
1048 : !-----------------------------------------------------------------------
1049 : subroutine outfld_mag( fldname, array )
1050 : use edyn_mpi, only: omlon1, mlon0, mlon1, mlat0, mlat1, mlev0, mlev1
1051 :
1052 : character(len=*), intent(in) :: fldname
1053 : real(r8), intent(in) :: array(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1)
1054 :
1055 : integer :: j
1056 :
1057 : do j = mlat0,mlat1
1058 : call outfld( fldname, array(mlon0:omlon1,j,mlev0:mlev1),omlon1-mlon0+1,j)
1059 : end do
1060 :
1061 : end subroutine outfld_mag
1062 : !-----------------------------------------------------------------------
1063 0 : subroutine outfld_mag2d( fldname, array )
1064 : use edyn_mpi, only: mlon0, mlon1, mlat0, mlat1
1065 :
1066 : character(len=*), intent(in) :: fldname
1067 : real(r8), intent(in) :: array(mlon0:mlon1,mlat0:mlat1)
1068 :
1069 : integer :: j
1070 :
1071 0 : do j = mlat0,mlat1
1072 0 : call outfld( fldname, array(mlon0:mlon1,j), mlon1-mlon0+1, j )
1073 : end do
1074 :
1075 0 : end subroutine outfld_mag2d
1076 :
1077 : end module dpie_coupling
|
