Line data Source code
1 : module edynamo
2 : !
3 : !-----------------------------------------------------------------------
4 : ! Purpose:
5 : ! Electro-dynamo module
6 : !-----------------------------------------------------------------------
7 : !
8 : use shr_kind_mod, only: r8 => shr_kind_r8 ! 8-byte reals
9 : use cam_logfile, only: iulog
10 : use cam_abortutils, only: endrun
11 : use spmd_utils, only: masterproc
12 : use edyn_params, only: finit ! initialization value
13 : use edyn_maggrid, only: nmlon, nmlonp1, nmlat, nmlath, nmlev
14 : use edyn_mpi, only: mlon0, mlon1, omlon1, mytid, mlat0, mlat1
15 : use edyn_mpi, only: mlev0, mlev1, lon0, lon1, lat0, lat1
16 : use edyn_solve, only: solve_edyn
17 : use time_manager, only: get_nstep ! for debug
18 : use cam_history, only: outfld, hist_fld_active
19 : use savefield_waccm, only: savefld_waccm
20 :
21 : implicit none
22 : save
23 : private
24 :
25 : integer :: nstep
26 : !
27 : ! 2d coefficients and RHS terms for PDE on magnetic subdomains
28 : ! (including halo points).
29 : ! If use_time3d_integ==.true., these will be input from time3d
30 : ! (see use-association in time3d.F90)
31 : !
32 : real(r8), allocatable, dimension(:,:) :: &
33 : zigm11, & ! sigma11*cos(theta0)
34 : zigmc, & ! sigmac
35 : zigm1, & ! for Hall conductance diagnostic
36 : zigm2, & ! sigma2
37 : zigm22, & ! sigma22/cos(theta0)
38 : rim1,rim2, & ! see description in comment below
39 : rhs, & ! right-hand side of PDE
40 : phimsolv, & ! solution direct from solver (nhem only)
41 : phim2d ! solution with phihm and both nhem and shem
42 : !
43 : ! 3d potential and electric field on mag subdomains (see sub pthreed):
44 : ! (mlon0:mlon1,mlat0:mlat1,mlev0:mlev1)
45 : ! Electric potential and field components are output fields of edynamo
46 : ! (later, these can be output arguments of the main driver, sub dynamo)
47 : !
48 : real(r8), allocatable, dimension(:,:,:) :: &
49 : phim3d, & ! 3d electric potential
50 : ed13d,ed23d, & ! 3d electric field for current calculations
51 : ephi3d, & ! 3d eastward electric field
52 : elam3d, & ! 3d equatorward electric field
53 : emz3d, & ! 3d upward electric field
54 : zpot_mag, &
55 : zpotm3d ! 3d geopotential (values at all levels)
56 : !
57 : ! 3d ion drift velocities on geographic grid (output):
58 : !
59 : ! real(r8), allocatable, dimension(:,:,:),save,target :: & ! (nlev,lon0:lon1,lat0:lat1)
60 : ! ui, & ! zonal ion drift
61 : ! vi, & ! meridional ion drift
62 : ! wi ! vertical ion drift
63 : !
64 : ! 3d electric field on geographic subdomains (see sub pefield):
65 : ! (nlev,lon0-2,lon1+2,lat0:lat1)
66 : real(r8), allocatable, dimension(:,:,:) :: ex,ey,ez
67 : !
68 : ! 3d electric potential on geographic subdomains (lon0:lon1,lat0:lat1,nlevp1)
69 : ! This will be regridded from phim3d for output to history files.
70 : real(r8), allocatable, dimension(:,:,:) :: phig3d ! (lon0:lon1,lat0:lat1,nlevp1)
71 : real(r8), allocatable, dimension(:,:,:) :: poten ! (nlevp1,lon0:lon1,lat0:lat1)
72 : !
73 : ! Fields at mag equator:
74 : !
75 : real(r8), allocatable, dimension(:,:) :: & ! (mlon0:mlon1,nmlev)
76 : ped_meq, hal_meq, adotv1_meq, adotv2_meq
77 : real(r8), allocatable, dimension(:,:,:) :: & ! (mlon0:mlon1,nmlev,4)
78 : fmeq_out
79 : real(r8), allocatable, dimension(:,:,:,:) :: & ! (mlon0:mlon1,mlat0:mlat1,nmlev,4)
80 : fmeq_in
81 : !
82 : ! Global longitude values near mag equator and poles for complete_integrals and rhs.
83 : ! These are declared in module data because they are used by subs complete_integrals
84 : ! and rhspde. The nf2d 7 fields are: zigm11,zigm22,zigmc,zigm1,zigm2,rim1,rim2,
85 : ! order is important (see feq_jpm1 and fpole_jpm2)!
86 : !
87 : integer, parameter :: nf2d=7 ! 7 2d fields
88 : real(r8), allocatable :: feq_jpm1(:,:,:) ! 7 fields at 2 lats (eq-1, eq+1)
89 : real(r8), allocatable :: fpole_jpm2(:,:,:) ! fields at S pole+1,2 and N pole-1,2
90 :
91 : real(r8), allocatable :: unitvm(:)
92 : !
93 : ! ed1,ed2: 2d electric field output on mag grid:
94 : ! (use-associated by dpie_coupling)
95 : !
96 : real(r8), allocatable, dimension(:,:) :: ed1, ed2 ! (mlon0-1:mlon1+1,mlat0-1:mlat1+1)
97 : !
98 : ! Global inputs to time3d: Note dimension order switch:
99 : ! edynamo has subdomains (mlon,mlat), whereas time3d has global (nmlat,nmlonp1)
100 : ! These are use-associated by time3d, and are init to zero in edyn_init.
101 : !
102 : real(r8), allocatable, dimension(:,:) :: ed1_glb, ed2_glb
103 : logical :: debug = .false. ! set true for prints to stdout at each call
104 :
105 : logical, public :: debug_hist = .false.
106 :
107 : public :: alloc_edyn, ed1, ed2, ed1_glb, ed2_glb
108 : public :: zigm11, zigmc, zigm2, zigm22, rim1, rim2
109 : public :: dynamo
110 :
111 : contains
112 : !-----------------------------------------------------------------------
113 7296 : subroutine dynamo( zpot_mag_in, ped_mag, hall_mag, adotv1_mag, adotv2_mag, adota1_mag, &
114 7296 : adota2_mag, a1dta2_mag,be3_mag, sini_mag, zpot, &
115 7296 : ui, vi, wi, lon0,lon1, lat0,lat1, lev0,lev1, do_integrals )
116 : use edyn_mpi, only: &
117 : mp_mag_halos, & ! set magnetic halo points
118 : mp_scatter_phim ! scatter solution to slave tasks
119 : use edyn_solve, only: rim_glb ! pde solver output (nmlonp1,nmlat,2)
120 : !
121 : ! Main driver for edynamo.
122 : ! Note alloc_edyn and esmf_init are called from edyn_init.
123 : !
124 : ! Args:
125 : integer,intent(in) :: & ! geographic subdomain
126 : lon0, lon1, & ! first,last longitude indices of geographic subdomain
127 : lat0, lat1, & ! first,last latitude indices of geographic subdomain
128 : lev0, lev1 ! first,last level indices (not distributed)
129 : !
130 : ! Inputs :
131 : !
132 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1), intent(in) :: &
133 : zpot_mag_in, & ! geopotential (cm)
134 : ped_mag, & ! pedersen conductivity (S/m)
135 : hall_mag, & ! hall conductivity (S/m)
136 : adotv1_mag,& ! ue1 (m/s)
137 : adotv2_mag ! ue2 (m/s)
138 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1), intent(in) :: &
139 : adota1_mag, &
140 : adota2_mag, &
141 : a1dta2_mag, &
142 : be3_mag, &
143 : sini_mag
144 :
145 : ! inputs on geographic (oplus) grid
146 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1), intent(in) :: &
147 : zpot ! geopotential (cm)
148 :
149 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1), intent(out) :: &
150 : ui, & ! zonal ion drift (cm/s)
151 : vi, & ! meridional ion drift (cm/s)
152 : wi ! vertical ion drift (cm/s)
153 :
154 : logical,intent(in) :: do_integrals
155 :
156 7296 : if (debug) then
157 0 : nstep = get_nstep()
158 0 : write(iulog,"(a,i5,a,l1)") 'Enter dynamo: nstep=', nstep, &
159 0 : ', do_integrals=', do_integrals
160 : end if
161 :
162 : !
163 : ! Regrid input fields from geographic to magnetic, and calculate
164 : ! some additional fields. If conductances are passed in from
165 : ! time3d (.not.do_integrals), then we do not need these inputs.
166 : !
167 7296 : if (do_integrals) then
168 7296 : call dynamo_set_data( zpot_mag_in, ped_mag, hall_mag, adotv1_mag, adotv2_mag )
169 7296 : if (debug) then
170 0 : write(iulog,"('edynamo debug: after dynamo_input')")
171 : end if
172 : end if
173 :
174 : !
175 : ! Fieldline integration:
176 : !
177 : ! If *not* doing fieldline integrations, then global conductances
178 : ! were passed in to the driver from time3d, and transformed from
179 : ! (nmlat,nmlonp1) to (nmlonp1,nmlat), defining zigmxx and rim1,2
180 : ! for the solver.
181 : !
182 7296 : if (do_integrals) then
183 : call fieldline_integrals(ped_mag, hall_mag, adotv1_mag, adotv2_mag, &
184 7296 : adota1_mag, adota2_mag, a1dta2_mag, be3_mag, sini_mag)
185 : end if
186 : !
187 : ! Equatorial and polar values, hemisphere folding:
188 : ! (these will be time3d integrations if do_integrals==.false.)
189 : !
190 7296 : call complete_integrals()
191 7296 : if (debug) then
192 0 : write(iulog,"('edynamo debug: after complete_integrals')")
193 : end if
194 : !
195 : ! Calculate right-hand side on mag subdomains:
196 : ! (mag halos are needed in rim1,2 for rhs calculation)
197 : !
198 7296 : call mp_mag_halos(rim1,mlon0,mlon1,mlat0,mlat1,1)
199 7296 : call mp_mag_halos(rim2,mlon0,mlon1,mlat0,mlat1,1)
200 7296 : call rhspde
201 7296 : if (debug) then
202 0 : write(iulog,"('edynamo debug: after rhspde')")
203 : end if
204 : !
205 : ! Gather needed arrays to root task for the serial solver:
206 : !
207 7296 : call gather_edyn()
208 7296 : if (debug) then
209 0 : write(iulog,"('edynamo debug: after gather_edyn')")
210 : end if
211 : !
212 : ! Root task now sets up stencils and calls the PDE solver:
213 : !
214 7296 : if (debug) then
215 0 : write(iulog,"('edynamo debug: call solve_edyn (master only)')")
216 : end if
217 7296 : if (mytid == 0) then
218 19 : call solve_edyn()
219 : end if
220 7296 : if (debug) write(iulog,"('edynamo debug: after solve_edyn (master only)')")
221 : !
222 : ! rim1 after solver is needed for highlat_poten. rim_glb is distributed
223 : ! to subdomains as rim1, and mag halos set. This will overwrite rim1 from
224 : ! fieldline_integrals, complete_integrals, etc.
225 : !
226 185991 : call mp_scatter_phim(rim_glb(:,:,1),rim1(mlon0:mlon1,mlat0:mlat1))
227 7296 : if (debug) then
228 0 : write(iulog,"('edynamo debug: after mp_scatter_phim')")
229 : end if
230 :
231 7296 : call mp_mag_halos(rim1,mlon0,mlon1,mlat0,mlat1,1)
232 7296 : if (debug) then
233 0 : write(iulog,"('edynamo debug: after mp_mag_halos')")
234 : end if
235 : !
236 : ! Add high latitude potential from empirical model (heelis or weimer)
237 : ! to solution rim1, defining phim2d on mag subdomains.
238 : !
239 7296 : call highlat_poten()
240 7296 : if (debug) then
241 0 : write(iulog,"('edynamo debug: after highlat_poten')")
242 : end if
243 : !
244 : ! Expand phim2d to phim3d, first setting mag halos in phim2d from
245 : ! hightlat_poten. phim3d will then be the final potential from pdynamo.
246 : !
247 7296 : call mp_mag_halos(phim2d,mlon0,mlon1,mlat0,mlat1,1)
248 :
249 7296 : call pthreed()
250 7296 : if (debug) then
251 0 : write(iulog,"('edynamo debug: after pthreed')")
252 : end if
253 : !
254 : ! Convert electric field to geographic grid:
255 7296 : call pefield()
256 7296 : if (debug) then
257 0 : write(iulog,"('edynamo debug: after pefield')")
258 : end if
259 :
260 : !
261 : ! Calculate ion drift velocities:
262 : !
263 :
264 7296 : call ionvel(zpot,ui,vi,wi, lon0,lon1, lat0,lat1, lev0,lev1)
265 7296 : if (debug) then
266 0 : write(iulog,"('edynamo debug: after ionvel')")
267 : end if
268 :
269 7296 : end subroutine dynamo
270 : !-----------------------------------------------------------------------
271 7296 : subroutine dynamo_set_data( zpot_mag_in, ped_mag, hall_mag, adotv1_mag, adotv2_mag )
272 : !
273 : !
274 : use edyn_params, only: h0, kbotdyn
275 : use edyn_mpi, only: mp_mageq ! get global values at mag equator
276 :
277 : !
278 : ! Args: Input fields on geographic grid:
279 : !
280 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),intent(in) :: &
281 : zpot_mag_in,&! cm
282 : ped_mag, & ! pedersen conductivity (S/m)
283 : hall_mag, & ! hall conductivity (S/m)
284 : adotv1_mag,& ! ue1 (m/s)
285 : adotv2_mag ! ue2 (m/s)
286 : !
287 : ! Local:
288 : !
289 : integer :: j, i, k
290 : !
291 7296 : if (debug .and. masterproc) then
292 0 : write(iulog,"('dynamo_input after savefld_waccm calls')")
293 0 : write(iulog,"('dynamo_input: kbotdyn=',i4)") kbotdyn
294 : end if
295 : !
296 : ! fmeq_in are input fields on 3d mag subdomains.
297 : ! allocate(fmeq_in(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1,4)
298 : !
299 23237646 : fmeq_in(:,:,:,1) = ped_mag(:,:,:)
300 23237646 : fmeq_in(:,:,:,2) = hall_mag(:,:,:)
301 23237646 : fmeq_in(:,:,:,3) = adotv1_mag(:,:,:)
302 23237646 : fmeq_in(:,:,:,4) = adotv2_mag(:,:,:)
303 : !
304 : ! Tasks w/ mag equator send eq data(i,k) to other tasks in their tidi:
305 : !
306 7296 : call mp_mageq(fmeq_in, fmeq_out, 4, mlon0, mlon1, mlat0, mlat1, nmlev)
307 : !
308 : ! Output arrays now have mag equator data on longitude subdomain
309 : ! and full column (mlon0:mlon1,nmlev)
310 : ! These will be used in fieldline_integrals.
311 : !
312 7358016 : ped_meq(:,:) = fmeq_out(:,:,1)
313 7358016 : hal_meq(:,:) = fmeq_out(:,:,2)
314 7358016 : adotv1_meq(:,:) = fmeq_out(:,:,3)
315 7358016 : adotv2_meq(:,:) = fmeq_out(:,:,4)
316 :
317 23237646 : zpot_mag(:,:,:) = zpot_mag_in(:,:,:)
318 : !
319 : ! Save geopotential on magnetic grid in zpotm3d, then
320 : ! limit max zpot_mag to h0 for use in fieldline integrals
321 : ! and pthreed. This should set zpot_mag to constant h0
322 : ! below kbotdyn. It is not necessary to set poles of zpotm3d
323 : ! since sub pthreed does not reference the poles of zpotm3d.
324 : !
325 955776 : do k = mlev0, mlev1
326 3830856 : do j = mlat0, mlat1
327 23230350 : do i=mlon0,mlon1
328 19406790 : zpotm3d(i,j,k) = zpot_mag(i,j,k)
329 22281870 : if (zpot_mag(i,j,k) < h0) then
330 11039926 : zpot_mag(i,j,k) = h0
331 : end if
332 : end do
333 : end do
334 : end do
335 29412 : do j = mlat0, mlat1
336 22064396 : call outfld('PED_MAG',ped_mag(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
337 22064396 : call outfld('HAL_MAG',hall_mag(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
338 22064396 : call outfld('ZPOT_MAG',zpot_mag(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
339 29412 : if (debug_hist) then
340 0 : call outfld('ADOTV1_MAG',adotv1_mag(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
341 0 : call outfld('ADOTV2_MAG',adotv2_mag(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
342 : endif
343 : end do
344 7296 : end subroutine dynamo_set_data
345 : !-----------------------------------------------------------------------
346 : !-----------------------------------------------------------------------
347 768 : subroutine alloc_edyn
348 : use edyn_geogrid, only: nlev
349 : !
350 : ! Allocate and initialize arrays for parallel dynamo (module data)
351 : ! (called once per run)
352 : !
353 : integer :: istat
354 : integer :: mlon00,mlon11,mlat00,mlat11
355 : !
356 768 : mlon00=mlon0-1 ; mlon11=mlon1+1
357 768 : mlat00=mlat0-1 ; mlat11=mlat1+1
358 : !
359 : ! 2d fields on mag subdomains (i,j):
360 : ! Certain fields are allocated with halos mlon0-1:mlon1+1,mlat0-1:mlat1+1
361 : !
362 3072 : allocate(zigm11(mlon00:mlon11,mlat00:mlat11),stat=istat)
363 768 : if (istat /= 0) call endrun('alloc_edyn: zigm11')
364 38442 : zigm11 = finit
365 2304 : allocate(zigmc(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
366 768 : if (istat /= 0) call endrun('alloc_edyn: zigmc')
367 38442 : zigmc = finit
368 2304 : allocate(zigm1(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
369 768 : if (istat /= 0) call endrun('alloc_edyn: zigm1')
370 38442 : zigm1 = finit
371 2304 : allocate(zigm2(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
372 768 : if (istat /= 0) call endrun('alloc_edyn: zigm2')
373 38442 : zigm2 = finit
374 2304 : allocate(zigm22(mlon00:mlon11,mlat00:mlat11),stat=istat)
375 768 : if (istat /= 0) call endrun('alloc_edyn: zigm22')
376 38442 : zigm22 = finit
377 2304 : allocate(rhs(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
378 768 : if (istat /= 0) call endrun('alloc_edyn: rhs')
379 38442 : rhs = finit
380 2304 : allocate(rim1(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
381 768 : if (istat /= 0) call endrun('alloc_edyn: rim1')
382 38442 : rim1 = finit
383 2304 : allocate(rim2(mlon00:mlon11,mlat00:mlat11) ,stat=istat)
384 768 : if (istat /= 0) call endrun('alloc_edyn: rim2')
385 38442 : rim2 = finit
386 2304 : allocate(phimsolv(mlon00:mlon11,mlat00:mlat11),stat=istat)
387 768 : if (istat /= 0) call endrun('alloc_edyn: phimsolv')
388 38442 : phimsolv = finit
389 2304 : allocate(phim2d(mlon00:mlon11,mlat00:mlat11),stat=istat)
390 768 : if (istat /= 0) call endrun('alloc_edyn: phim2d')
391 38442 : phim2d = finit
392 : !
393 : ! 3d phim and electric field on mag subdomains:
394 3840 : allocate(phim3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
395 768 : if (istat /= 0) call endrun('alloc_edyn: phim3d')
396 2446068 : phim3d = finit
397 3840 : allocate(ed13d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
398 768 : if (istat /= 0) call endrun('alloc_edyn: ed13d')
399 2446068 : ed13d = finit
400 3840 : allocate(ed23d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
401 768 : if (istat /= 0) call endrun('alloc_edyn: ed23d')
402 2446068 : ed23d = finit
403 3840 : allocate(ephi3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
404 768 : if (istat /= 0) call endrun('alloc_edyn: ephi3d')
405 2446068 : ephi3d = finit
406 3840 : allocate(elam3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
407 768 : if (istat /= 0) call endrun('alloc_edyn: elam3d')
408 2446068 : elam3d = finit
409 3840 : allocate(emz3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
410 768 : if (istat /= 0) call endrun('alloc_edyn: emz3d')
411 2446068 : emz3d = finit
412 3840 : allocate(zpotm3d(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
413 768 : if (istat /= 0) call endrun('alloc_edyn: zpotm3d')
414 2446068 : zpotm3d = finit
415 3840 : allocate(zpot_mag(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1),stat=istat)
416 768 : if (istat /= 0) call endrun('alloc_edyn: zpot_mag')
417 2446068 : zpot_mag = finit
418 : !
419 : ! Fields at mag equator (subdomain longitudes and full column):
420 : !
421 3072 : allocate(ped_meq(mlon0:mlon1,mlev0:mlev1),stat=istat)
422 768 : if (istat /= 0) call endrun('alloc_edyn: ped_meq')
423 774528 : ped_meq = finit
424 3072 : allocate(hal_meq(mlon0:mlon1,mlev0:mlev1),stat=istat)
425 768 : if (istat /= 0) call endrun('alloc_edyn: hal_meq')
426 774528 : hal_meq = finit
427 3072 : allocate(adotv1_meq(mlon0:mlon1,mlev0:mlev1),stat=istat)
428 768 : if (istat /= 0) call endrun('alloc_edyn: adotv1_meq')
429 774528 : adotv1_meq = finit
430 3072 : allocate(adotv2_meq(mlon0:mlon1,mlev0:mlev1),stat=istat)
431 768 : if (istat /= 0) call endrun('alloc_edyn: adotv2_meq')
432 774528 : adotv2_meq = finit
433 : !
434 : ! Fields input to mp_mageq (4 fields at full mag subdomain i,j,k):
435 : !
436 4608 : allocate(fmeq_in(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1,4),stat=istat)
437 768 : if (istat /= 0) call endrun('alloc_edyn: fmeq_in')
438 9785040 : fmeq_in = finit
439 : !
440 : ! Fields output by mp_mageq (4 fields at mag subdomain i,k)
441 : !
442 3840 : allocate(fmeq_out(mlon0:mlon1,mlev0:mlev1,4),stat=istat)
443 768 : if (istat /= 0) call endrun('alloc_edyn: fmeq_out')
444 3098880 : fmeq_out = finit
445 : !
446 : ! 3d electric field on geographic subdomains, with halos:
447 : !
448 3840 : allocate(ex(nlev,lon0:lon1,lat0:lat1),stat=istat)
449 768 : if (istat /= 0) call endrun('alloc_edyn: ex')
450 3840 : allocate(ey(nlev,lon0:lon1,lat0:lat1),stat=istat)
451 768 : if (istat /= 0) call endrun('alloc_edyn: ey')
452 3840 : allocate(ez(nlev,lon0:lon1,lat0:lat1),stat=istat)
453 768 : if (istat /= 0) call endrun('alloc_edyn: ez')
454 10873344 : ex=finit ; ey=finit ; ez=finit
455 : !
456 : ! 3d electric potential on geographic subdomains (k,i,j):
457 3840 : allocate(poten(nlev,lon0:lon1,lat0:lat1),stat=istat)
458 768 : if (istat /= 0) call endrun('alloc_pdyn: poten')
459 3624960 : poten = finit
460 : !
461 : ! 3d electric potential on geographic subdomains ((i,j,k) regridded from phim3d):
462 3840 : allocate(phig3d(lon0:lon1,lat0:lat1,mlev0:mlev1),stat=istat)
463 768 : if (istat /= 0) call endrun('alloc_pdyn: phig3d')
464 3994368 : phig3d = finit
465 : !
466 : ! 2d electric field components on mag grid (these may be input to time3d):
467 : ! real(r8), dimension(:,:) :: ed1,ed2 ! (mlon0-1:mlon1+1,mlat0-1:mlat1+1)
468 : !
469 3072 : allocate(ed1(mlon0-1:mlon1+1,mlat0-1:mlat1+1),stat=istat)
470 768 : if (istat /= 0) call endrun('alloc_edyn: ed1')
471 38442 : ed1 = finit
472 :
473 3072 : allocate(ed2(mlon0-1:mlon1+1,mlat0-1:mlat1+1),stat=istat)
474 768 : if (istat /= 0) call endrun('alloc_edyn: ed2')
475 38442 : ed2 = finit
476 :
477 2304 : allocate(unitvm(nmlon))
478 62208 : unitvm = 1._r8
479 :
480 3072 : allocate(feq_jpm1(nmlonp1,2,nf2d))
481 3072 : allocate(fpole_jpm2(nmlonp1,4,nf2d))
482 4608 : allocate(ed1_glb(nmlat,nmlonp1), ed2_glb(nmlat,nmlonp1))
483 :
484 768 : end subroutine alloc_edyn
485 : !-----------------------------------------------------------------------
486 7296 : subroutine fieldline_integrals( ped_mag, hal_mag, adotv1_mag, adotv2_mag, &
487 7296 : adota1_mag, adota2_mag, a1dta2_mag, be3_mag, sini_mag )
488 : !
489 : ! Integrate along magnetic field lines, saving conductances and rims.
490 : !
491 : use edyn_params, only: r0,h0,finit,kbotdyn
492 : use edyn_maggrid, only: ylatm
493 : !
494 : ! Args:
495 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1,mlev0:mlev1), intent(in) :: &
496 : ped_mag, & ! pedersen conductivity (S/m)
497 : hal_mag, & ! hall conductivity (S/m)
498 : adotv1_mag,& ! ue1 (m/s)
499 : adotv2_mag ! ue2 (m/s)
500 :
501 : real(r8), dimension(mlon0:mlon1,mlat0:mlat1), intent(in) :: &
502 : adota1_mag, &
503 : adota2_mag, &
504 : a1dta2_mag, &
505 : be3_mag, &
506 : sini_mag
507 :
508 : !
509 : ! Local:
510 : integer :: i,j,k
511 : real(r8) :: &
512 : sinlm, & ! sin(lam_m)
513 : clm2, & ! cos^2(lam_m)
514 : ra, & ! (R_E + H_0)/cos^2(lam_m)
515 : sqomrra, & ! sqrt(1/ R_0/R_A) = sin(lam_m)
516 : sqrra, & ! sqrt(R_0/R_A)
517 : afac, & ! 2*sqrt(R_A-R_0) (afac is NaN at mag poles)
518 : htfac ! sqrt(R_A -3/4*R_0)
519 :
520 : real(r8) :: rora,del,omdel,sig1,sig2,ue1,ue2
521 14592 : real(r8), dimension(mlon0:mlon1) :: aam
522 14592 : real(r8), dimension(mlon0:mlon1,mlev0:mlev1) :: rrm, &
523 14592 : rtramrm, htfunc, htfunc2
524 :
525 : !
526 : ! Initialize coefficients:
527 : !
528 365199 : zigm11 = finit
529 365199 : zigm22 = finit
530 365199 : zigm1 = finit
531 365199 : zigm2 = finit
532 365199 : zigmc = finit
533 365199 : rim1 = finit
534 365199 : rim2 = finit
535 : !
536 : ! Subdomain latitude scan for field line integrations:
537 : !
538 29412 : do j=mlat0,mlat1
539 : !
540 : ! Skip poles and equator:
541 22116 : if (j==1.or.j==nmlat.or.j==nmlath) cycle
542 :
543 21432 : sinlm = sin(ylatm(j))
544 21432 : clm2 = 1._r8 - sinlm*sinlm
545 21432 : ra = r0/clm2
546 21432 : sqomrra = sqrt(1._r8-r0/ra)
547 21432 : sqrra = sqrt(r0/ra)
548 21432 : afac = 2._r8*sqrt(ra-r0)
549 21432 : htfac = sqrt(ra-0.75_r8*r0)
550 166098 : do i=mlon0,mlon1
551 166098 : aam(i) = afac/abs(sini_mag(i,j))
552 : enddo
553 : !
554 : ! 2*sqrt( h_A - h_0 )/ |sin I_m | w.r to reference height A(h_R)
555 1500240 : do k=kbotdyn,mlev1
556 11482194 : do i=mlon0,mlon1
557 : !
558 : ! rr = r0+z-z0 radius of magnetic point
559 : ! (Note zpot_mag min was set to h0 in dynamo_inputs)
560 9981954 : rrm(i,k) = r0+zpot_mag(i,j,k)-h0
561 : !
562 : ! rtramr = ra-r if +ive, zero otherwise
563 9981954 : rtramrm(i,k) = max(0._r8,ra-rrm(i,k))
564 11460762 : rtramrm(i,k) = sqrt(rtramrm(i,k))
565 : enddo ! i=mlon0,mlon1
566 : enddo ! k=kbotdyn,mlev1
567 : !
568 : ! Interpolate to midpoints:
569 : ! htfunc = factor by which to multiply AAM(I)*d(sqrt(ra-r)) = ds
570 1478808 : do k=kbotdyn,mlev1-1
571 11316096 : do i=mlon0,mlon1
572 9837288 : rrm(i,k) = 0.5_r8*(rrm(i,k)+rrm(i,k+1))
573 9837288 : rtramrm(i,k) = rtramrm(i,k)-rtramrm(i,k+1)
574 9837288 : htfunc(i,k) = sqrt(ra-0.75_r8*rrm(i,k))/htfac
575 11294664 : htfunc2(i,k) = htfunc(i,k)**2
576 : enddo ! i=mlon0,mlon1
577 : enddo ! k=kbotdyn,mlev1
578 : !
579 : ! Compute integrals:
580 1478808 : do k=kbotdyn,mlev1-1
581 11316096 : do i=mlon0,mlon1
582 : !
583 : ! (R_E+h)/(R_E+h_A) < 1 -> h_A > h
584 9837288 : rora = min(1._r8,rrm(i,k)/ra)
585 : !
586 : ! (lam_m - lam) / lam_m =
587 : ! sqrt(1-r_0/r_A)sqrt(r/r_A) - sqrt(r_0/r_A)sqrt(1-r/r_A)
588 9837288 : del = (sqomrra*sqrt(rora)-sqrra*sqrt(1._r8-rora))/abs(ylatm(j))
589 9837288 : omdel = 1._r8 - del
590 : !
591 : ! Interpolate conductivities and winds in latitude along field line, assuming
592 : ! linear variation between foot of field line and magnetic equator.
593 : ! (For field lines other than those near the magnetic equator, del is nearly
594 : ! zero, so that the interpolated values are essentially the values for the
595 : ! latitude of the foot of the field line; inaccuracy of the assumption of
596 : ! linear variation is thus unimportant for these field lines.)
597 : !
598 : ! Here, mag equator ped_meq, etc. are from mp_mageq, called from dynamo inputs:
599 9837288 : sig1 = omdel*ped_mag(i,j,k) + del*ped_meq(i,k)
600 9837288 : sig2 = omdel*hal_mag(i,j,k) + del*hal_meq(i,k)
601 9837288 : ue1 = omdel*adotv1_mag(i,j,k) + del*adotv1_meq(i,k)
602 9837288 : ue2 = omdel*adotv2_mag(i,j,k) + del*adotv2_meq(i,k)
603 : !
604 : ! height varying factors: ds = aam*htfunc
605 : ! d_1^2/D = 1/htfunc * adota1_mag(i,j)
606 : ! d_2^2/D = htfunc * adota2_mag(i,j)
607 : ! d_1*d_2/D = 1 * a1dta2_mag(i,j)
608 : !
609 : ! zigm11: int (sigma_p*d_1^2/D) ds : d_1^2/D
610 : ! zigm22: int (sigma_p*d_2^2/D) ds : d_2^2/D
611 : !
612 9837288 : zigm11(i,j) = zigm11(i,j) + sig1*rtramrm(i,k)
613 9837288 : zigm22(i,j) = zigm22(i,j) + sig1*rtramrm(i,k)*htfunc2(i,k)
614 :
615 : !
616 : ! zigmc: int (sigma_p*d_1*d_2/D) ds
617 : ! zigm2: int (sigma_h) ds
618 : !
619 9837288 : zigmc(i,j) = zigmc(i,j) + sig1*rtramrm(i,k)*htfunc(i,k)
620 9837288 : zigm2(i,j) = zigm2(i,j) + sig2*rtramrm(i,k)*htfunc(i,k)
621 :
622 : ! zigm1: int(sigma_p) ds
623 9837288 : zigm1(i,j) = zigm1(i,j) + sig1*rtramrm(i,k)*htfunc(i,k)
624 : !
625 : ! rim1: int [sigma_p*d_1^2/D u_e2+(sigma_h-(sigma_p*d_1*d_2)/D) u_e1] ds
626 : ! rim2: int [(sigma_h+sigma_p*d_1*d_2/D) u_e2-sigma_p*d_2^2/D u_e1 ] ds
627 : !
628 0 : rim1(i,j) = rim1(i,j) + (sig1*adota1_mag(i,j)*ue2 + &
629 9837288 : (sig2 - sig1*a1dta2_mag(i,j))*htfunc(i,k)*ue1)*rtramrm(i,k)
630 0 : rim2(i,j) = rim2(i,j) + (sig1*adota2_mag(i,j)*htfunc2(i,k)* &
631 11294664 : ue1 - (sig2 + sig1*a1dta2_mag(i,j))*htfunc(i,k)*ue2)*rtramrm(i,k)
632 : enddo ! i=mlon0,mlon1
633 : enddo ! k=kbotdyn,mlev1
634 : !
635 : ! Complete calculation and place result in /coefm/ zigm's
636 : ! rim's are in A/m multiply by 1/100 to convert from [cm] to [m]
637 : !
638 : ! At this point,
639 : ! zigm11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im)
640 : ! zigm22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im)
641 : ! zigmc is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c
642 : ! zigm2 is int[sigma_h] ds, i.e. Sigma_h
643 : !
644 : ! rim1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)*
645 : ! B_e3 ds, i.e. K_(m phi)^D/abs(sin Im)
646 : ! rim2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)*
647 : ! B_e3 ds, K_(m lam)^D ( minus in northern hemisphere
648 : ! Change sign of RIM(2) in S. hemisphere to be compatible with transf
649 : ! At this point, rim2 is +-K_(m lam)^D
650 : !
651 173394 : do i = mlon0,mlon1
652 144666 : zigm11(i,j) = 1.e-2_r8*zigm11(i,j)*aam(i)*adota1_mag(i,j)
653 144666 : zigm22(i,j) = 1.e-2_r8*zigm22(i,j)*aam(i)*adota2_mag(i,j)
654 144666 : zigmc(i,j) = 1.e-2_r8*zigmc (i,j)*aam(i)*a1dta2_mag(i,j)
655 144666 : zigm2(i,j) = 1.e-2_r8*zigm2 (i,j)*aam(i)
656 144666 : zigm1(i,j) = 1.e-2_r8*zigm1 (i,j)*aam(i)
657 144666 : rim1(i,j) = 1.e-2_r8*rim1(i,j)*aam(i)*be3_mag(i,j)
658 166782 : rim2(i,j) = 1.e-2_r8*rim2(i,j)*aam(i)*be3_mag(i,j)
659 : enddo ! i = 1,nmlon
660 : enddo ! j=mlat0,mlat1 (without poles)
661 :
662 7296 : if (debug_hist) then
663 0 : call savefld_waccm(adota1_mag(mlon0:mlon1,mlat0:mlat1) ,'adota1_mag_a' ,1,mlon0,mlon1,mlat0,mlat1)
664 0 : call savefld_waccm(zigm11(mlon0:mlon1,mlat0:mlat1) ,'ZIGM11_a' ,1,mlon0,mlon1,mlat0,mlat1)
665 : endif
666 :
667 7296 : end subroutine fieldline_integrals
668 : !-----------------------------------------------------------------------
669 7296 : subroutine complete_integrals
670 : use edyn_mpi, only: mlat0, mlat1, mlon0, mlon1, mp_mageq_jpm1
671 : use edyn_mpi, only: mp_magpole_2d, mp_mag_foldhem, mp_mag_periodic_f2d
672 : use edyn_maggrid, only: rcos0s,dt1dts
673 : !
674 : ! Field line integrals for each hemisphere have been calculated in
675 : ! mag subdomains. Now, complete these arrays with equator and polar
676 : ! values, and sum the integrals from the 2 hemispheres.
677 : ! This is done by obtaining global mag fields via mpi exchange
678 : ! of mag subdomains, completing the global fields, and updating
679 : ! the subdomains.
680 : ! The 6 2d fields are: zigm11,zigm22,zigmc,zigm2,rim1,rim2
681 : !
682 : ! Local:
683 : integer :: i,j,ii,lonend
684 14592 : real(r8) :: fmsub(mlon0:mlon1,mlat0:mlat1,nf2d)
685 : real(r8) :: corfac
686 : real(r8) :: r8_nmlon
687 7296 : r8_nmlon = real(nmlon, r8)
688 : !
689 : ! For equatorial values, we need latitudes eq+1 and eq-1:
690 : ! Local feq_jpm1(nmlonp1,2,6) is returned by mp_mageq_jpm1,
691 : ! where the 2 dim contains lats nmlath-1, nmlath+1. These
692 : ! are global in lon, even tho each subd uses only its own i's.
693 : ! These mag equator values do not show up on plots because
694 : ! of the small factor .06 and .125.
695 : ! The 6 2d fields are: zigm11,zigm22,zigmc,zigm2,rim1,rim2
696 : ! Must specify mlon0:mlon1,mlat0:mlat1 because these fields
697 : ! were allocated to include single-point halos (these calls
698 : ! exclude the halo points):
699 : !
700 178695 : fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1)
701 178695 : fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1)
702 178695 : fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1)
703 178695 : fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1)
704 178695 : fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1)
705 178695 : fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1)
706 178695 : fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1)
707 :
708 7296 : call mp_mageq_jpm1(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1,feq_jpm1,nf2d)
709 : !
710 : ! From sub fieldline_integrals:
711 : ! zigm11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im)
712 : ! zigm22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im)
713 : ! zigmc is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c
714 : ! zigm2 is int[sigma_h] ds, i.e. Sigma_h
715 : !
716 : ! rim1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)*
717 : ! B_e3 ds, i.e. K_(m phi)^D/abs(sin Im)
718 : ! rim2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)*
719 : ! B_e3 ds, K_(m lam)^D ( minus in northern hemisphere
720 : ! Change sign of RIM(2) in S. hemisphere to be compatible with transf
721 : ! At this point, rim2 is +-K_(m lam)^D
722 : !
723 : ! Equatorial values:
724 : ! Assume that quantities primarily dependent on Pedersen conductivity
725 : ! have field-line integrals 1/4 as large as the averages for next-higher
726 : ! field lines; quantities primarily dependent on Hall conductivity
727 : ! have field-line integrals 0.12 as large as the averages for next-higher
728 : ! field lines. Exact values chosen should not be important for potential
729 : ! calculation, as long as they are physically reasonable and not too
730 : ! different from adjacent values.
731 : !
732 29412 : do j=mlat0,mlat1
733 29412 : if (j==nmlath) then ! mag equator
734 :
735 0 : zigm11(mlon0:mlon1,j) = .125_r8*(feq_jpm1(mlon0:mlon1,1,1)+ &
736 1767 : feq_jpm1(mlon0:mlon1,2,1))
737 0 : zigm22(mlon0:mlon1,j) = .125_r8*(feq_jpm1(mlon0:mlon1,1,2)+ &
738 1767 : feq_jpm1(mlon0:mlon1,2,2))
739 0 : zigmc (mlon0:mlon1,j) = .125_r8*(feq_jpm1(mlon0:mlon1,1,3)+ &
740 1767 : feq_jpm1(mlon0:mlon1,2,3))
741 0 : zigm2 (mlon0:mlon1,j) = .060_r8*(feq_jpm1(mlon0:mlon1,1,4)+ &
742 1767 : feq_jpm1(mlon0:mlon1,2,4))
743 0 : rim1 (mlon0:mlon1,j) = .060_r8*(feq_jpm1(mlon0:mlon1,1,5)+ &
744 1767 : feq_jpm1(mlon0:mlon1,2,5))
745 0 : rim2 (mlon0:mlon1,j) = .060_r8*(feq_jpm1(mlon0:mlon1,1,6)+ &
746 1767 : feq_jpm1(mlon0:mlon1,2,6))
747 0 : zigm1 (mlon0:mlon1,j) = .060_r8*(feq_jpm1(mlon0:mlon1,1,7)+ &
748 1767 : feq_jpm1(mlon0:mlon1,2,7))
749 : !
750 : ! Include the boundary condition at the equator eq.(5.30) in
751 : ! Richmond (1995) Ionospheric Electrodynamics use. Mag. Apex Coord.
752 : ! J.Geomag.Geoelectr. 47,191-212
753 : ! Sig_phiphi/abs(sin Im) = 0.5*Sig_cowling/abs(sin Im)
754 : ! = 0.5/abs(sin Im)*(Sig_phiphi - Sig_philam*sig_lamphi/Sig_lamlam)
755 : ! = 0.5/abs(sin Im)*(Sig_phiphi + (Sig_h-sig_c)*(Sig_h+sig_c)/Sig_lamlam)
756 : ! rim(1) / |sin I_m| = I_1 = R/2*(K_mphi - Sig_philam/Sig_lamlam*K_mlam)
757 : !
758 1767 : do i=mlon0,mlon1
759 0 : zigm11(i,j) = zigm11(i,j)+(zigm2(i,j)-zigmc(i,j))* &
760 1539 : (zigm2(i,j)+zigmc(i,j))/zigm22(i,j)
761 0 : rim1(i,j) = rim1(i,j) - (zigm2(i,j)-zigmc(i,j))/ &
762 1767 : zigm22(i,j)*rim2(i,j)
763 : enddo ! i=mlon0,mlon1
764 : endif ! j at equator
765 : enddo ! j=mlat0,mlat1
766 : !
767 : ! Using notation of Richmond (1995) on right-hand side below:
768 : ! Sigma_(phi phi) = zigm11*abs(sin I_m)
769 : ! Sigma_(lam lam) = zigm22/abs(sin I_m)
770 : ! Sigma_(phi lam) = +-(zigm2-zigmc)
771 : ! Sigma_(lam phi) = -+(zigm2+zigmc)
772 : ! K_(m phi)^D = rim(1)*abs(sin I_m)
773 : ! K_(m lam)^D = +-rim(2)
774 : !
775 : ! Transforming PDE from original apex (theta_a) to new apex grid (theta_0)
776 : ! which is equally spaced in magnetic latitude
777 : ! SCALE quantities to modified (0) magnetic latitude system, multiplying or
778 : ! dividing by abs(sin I_m) [inverse contained in DT1DTS] as necessary.
779 : ! Sign of K_(m lam)^D in southern hemisphere remains reversed.
780 : ! for the mixed terms the transformation from the integration and differentiation
781 : ! canceled out (zigmc, zigm2)
782 : ! DT1DTS : d theta_0/ d theta_a / abs(sin I_m)
783 : ! RCOS0S : cos(theta_0)/ cos(theta_a)
784 : !
785 : ! corfac: abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a)
786 : ! zigm11: abs(I_m)*d theta_a/d theta_0 * cos(theta_0)/ cos(theta_a)
787 : ! zigm22: 1/abs(I_m)*d theta_0/d theta_a * cos(theta_a)/ cos(theta_0)
788 : ! rim(1): abs(I_m)*d theta_a/d theta_0
789 : ! rim(2): cos(theta_a)/ cos(theta_0)
790 : !
791 29412 : do j=mlat0,mlat1
792 22116 : if (j==1.or.j==nmlat) cycle ! skip poles
793 21660 : corfac = rcos0s(j)/dt1dts(j)
794 175161 : do i=mlon0,mlon1
795 146205 : zigm11(i,j) = zigm11(i,j)*corfac
796 146205 : zigm22(i,j) = zigm22(i,j)/corfac
797 146205 : rim1(i,j) = rim1(i,j)/dt1dts(j)
798 168321 : rim2(i,j) = rim2(i,j)/rcos0s(j)
799 : enddo
800 : enddo
801 : !
802 : ! For polar values, we need south pole plus 1 and 2 (j==2,3),
803 : ! and north pole minus 1 and 2 (j==nmlat-1,nmlat-2). These
804 : ! are returned by sub mp_magpole_jpm2 (mpi.F):
805 : ! Must specify (mlon0:mlon1,mlat0:mlat1) because zigmxx and rims
806 : ! are allocated to include halo cells.
807 : !
808 178695 : fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1)
809 185991 : fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1)
810 185991 : fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1)
811 185991 : fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1)
812 185991 : fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1)
813 185991 : fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1)
814 185991 : fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1)
815 : !
816 : ! mp_magpole_2d returns fpole_jpm2(nmlonp1,1->4,nf) as:
817 : ! 1: j = 2 (spole+1)
818 : ! 2: j = 3 (spole+2)
819 : ! 3: j = nmlat-1 (npole-1)
820 : ! 4: j = nmlat-2 (npole-2)
821 : !
822 : call mp_magpole_2d(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1, &
823 7296 : 1,nmlat,fpole_jpm2,nf2d)
824 : !
825 : ! the PDE is divided by 1/ DT0DTS
826 : ! Sigma_(phi phi) = zigm11/ rcos0s * dt0dts
827 : ! Sigma_(lam lam) = zigm22 * rcos0s / dt0dts
828 : ! Sigma_(phi lam) = +-(zigm2-zigmc)
829 : ! Sigma_(lam phi) = -+(zigm2+zigmc)
830 : ! K_(m phi)^D = rim(1) * dt0dts
831 : ! K_(m lam)^D = +-rim(2) * rcos0s
832 : !
833 : ! Compute polar values for the conductances, 4th order interpolation:
834 : !
835 29412 : do j=mlat0,mlat1
836 : !
837 : ! South pole:
838 29412 : if (j==1) then ! south pole (use fpole_jpm2(nmlon,1->2,nf)
839 0 : zigm11(mlon0,j)=(4._r8* &
840 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,1,1))- &
841 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,2,1)))/ &
842 37164 : (3._r8*r8_nmlon)
843 0 : zigm22(mlon0,j)=(4._r8* &
844 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,1,2))- &
845 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,2,2)))/ &
846 37164 : (3._r8*r8_nmlon)
847 0 : zigmc (mlon0,j)=(4._r8* &
848 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,1,3))- &
849 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,2,3)))/ &
850 37164 : (3._r8*r8_nmlon)
851 0 : zigm2 (mlon0,j)=(4._r8* &
852 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,1,4))- &
853 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,2,4)))/ &
854 37164 : (3._r8*r8_nmlon)
855 0 : zigm1(mlon0,j) = (4._r8* &
856 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,1,7))- &
857 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,2,7)))/ &
858 37164 : (3._r8*r8_nmlon)
859 : !
860 : ! Extend south pole over longitude:
861 1539 : do i=mlon0+1,mlon1
862 1311 : zigm11(i,j) = zigm11(mlon0,j)
863 1311 : zigm22(i,j) = zigm22(mlon0,j)
864 1311 : zigmc (i,j) = zigmc (mlon0,j)
865 1311 : zigm2 (i,j) = zigm2 (mlon0,j)
866 1539 : zigm1 (i,j) = zigm1 (mlon0,j)
867 : enddo ! i=mlon0,mlon1
868 : !
869 : ! RHS vector (I_1,I_2): average over south pole:
870 : ! (use fpole_jpm2(i,1,nf), i.e. j==2, and lons across the pole)
871 228 : lonend = mlon1
872 228 : if (mlon1==nmlonp1) lonend = mlon1-1
873 1748 : do i=mlon0,lonend
874 1520 : ii = 1+mod(i-1+nmlon/2,nmlon)
875 1520 : rim1(i,j) = 0.5_r8*(fpole_jpm2(i,1,5)-fpole_jpm2(ii,1,5))
876 1748 : rim2(i,j) = 0.5_r8*(fpole_jpm2(i,1,6)-fpole_jpm2(ii,1,6))
877 : enddo
878 : !
879 : ! North pole:
880 21888 : elseif (j==nmlat) then ! north pole (use fpole_jpm2(nmlon,3->4,1,nf)
881 0 : zigm11(mlon0,j)=(4._r8* &
882 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,3,1))- &
883 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,4,1)))/ &
884 37164 : (3._r8*r8_nmlon)
885 0 : zigm22(mlon0,j)=(4._r8* &
886 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,3,2))- &
887 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,4,2)))/ &
888 37164 : (3._r8*r8_nmlon)
889 0 : zigmc (mlon0,j)=(4._r8* &
890 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,3,3))- &
891 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,4,3)))/ &
892 37164 : (3._r8*r8_nmlon)
893 0 : zigm2 (mlon0,j)=(4._r8* &
894 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,3,4))- &
895 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,4,4)))/ &
896 37164 : (3._r8*r8_nmlon)
897 0 : zigm1(mlon0,j) = (4._r8* &
898 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,3,7))- &
899 228 : dot_product(unitvm,fpole_jpm2(1:nmlon,4,7)))/ &
900 37164 : (3._r8*r8_nmlon)
901 : !
902 : ! Extend north pole over longitude:
903 1539 : do i=mlon0+1,mlon1
904 1311 : zigm11(i,j) = zigm11(mlon0,j)
905 1311 : zigm22(i,j) = zigm22(mlon0,j)
906 1311 : zigmc (i,j) = zigmc (mlon0,j)
907 1311 : zigm2 (i,j) = zigm2 (mlon0,j)
908 1539 : zigm1 (i,j) = zigm1 (mlon0,j)
909 : enddo ! i=mlon0,mlon1
910 : !
911 : ! RHS vector (I_1,I_2): average over north pole:
912 : ! (use fpole_jpm2(i,3,nf), i.e. j==nmlat-1, and lons across the pole)
913 228 : lonend = mlon1
914 228 : if (mlon1==nmlonp1) lonend = mlon1-1
915 1748 : do i=mlon0,lonend
916 1520 : ii = 1+mod(i-1+nmlon/2,nmlon)
917 1520 : rim1(i,j) = 0.5_r8*(fpole_jpm2(i,3,5)-fpole_jpm2(ii,3,5))
918 1748 : rim2(i,j) = 0.5_r8*(fpole_jpm2(i,3,6)-fpole_jpm2(ii,3,6))
919 : enddo
920 : endif ! south or north pole
921 : enddo ! j=mlat0,mlat1
922 : !
923 : ! Fold south hemisphere over onto north, and set periodic points:
924 178695 : fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1)
925 185991 : fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1)
926 185991 : fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1)
927 185991 : fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1)
928 185991 : fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1)
929 185991 : fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1)
930 185991 : fmsub(:,:,7) = zigm1 (mlon0:mlon1,mlat0:mlat1)
931 :
932 7296 : call mp_mag_foldhem(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d)
933 7296 : call mp_mag_periodic_f2d(fmsub,mlon0,mlon1,mlat0,mlat1,nf2d)
934 :
935 178695 : zigm11(mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,1)
936 185991 : zigm22(mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,2)
937 185991 : zigmc (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,3)
938 185991 : zigm2 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,4)
939 185991 : rim1 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,5)
940 185991 : rim2 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,6)
941 185991 : zigm1 (mlon0:mlon1,mlat0:mlat1) = fmsub(:,:,7)
942 : !
943 : ! Reverse sign of zigmc in northern hemisphere.
944 29412 : do j=mlat0,mlat1
945 22116 : if (j >= nmlath) then
946 86583 : zigmc(mlon0:mlon1,j) = -zigmc(mlon0:mlon1,j)
947 : endif
948 22116 : if (debug_hist) then
949 0 : call outfld('EDYN_RIM1',rim1(mlon0:omlon1,j),omlon1-mlon0+1,j)
950 0 : call outfld('EDYN_RIM2',rim2(mlon0:omlon1,j),omlon1-mlon0+1,j)
951 : endif
952 22116 : call outfld('PED_CONDUCTANCE', zigm2(mlon0:omlon1,j),omlon1-mlon0+1,j)
953 29412 : call outfld('HALL_CONDUCTANCE',zigm1(mlon0:omlon1,j),omlon1-mlon0+1,j)
954 : enddo
955 :
956 7296 : if (debug.and.masterproc) then
957 0 : write(iulog,"('complete_integrals: nstep=',i4)") nstep
958 : write(iulog,"(' zigm11 min,max=',2e12.4)") &
959 0 : minval(zigm11(mlon0:mlon1,mlat0:mlat1)),maxval(zigm11(mlon0:mlon1,mlat0:mlat1))
960 : write(iulog,"(' zigm22 min,max=',2e12.4)") &
961 0 : minval(zigm22(mlon0:mlon1,mlat0:mlat1)),maxval(zigm22(mlon0:mlon1,mlat0:mlat1))
962 : write(iulog,"(' zigmc min,max=',2e12.4)") &
963 0 : minval(zigmc (mlon0:mlon1,mlat0:mlat1)),maxval(zigmc (mlon0:mlon1,mlat0:mlat1))
964 : write(iulog,"(' zigm2 min,max=',2e12.4)") &
965 0 : minval(zigm2 (mlon0:mlon1,mlat0:mlat1)),maxval(zigm2 (mlon0:mlon1,mlat0:mlat1))
966 : write(iulog,"(' rim1 min,max=',2e12.4)") &
967 0 : minval(rim1 (mlon0:mlon1,mlat0:mlat1)),maxval(rim1 (mlon0:mlon1,mlat0:mlat1))
968 : write(iulog,"(' rim2 min,max=',2e12.4)") &
969 0 : minval(rim2 (mlon0:mlon1,mlat0:mlat1)),maxval(rim2 (mlon0:mlon1,mlat0:mlat1))
970 : endif
971 :
972 7296 : if (debug_hist) then
973 0 : call savefld_waccm(zigm11(mlon0:mlon1,mlat0:mlat1) ,'EDYN_ZIGM11' ,1,mlon0,mlon1,mlat0,mlat1)
974 0 : call savefld_waccm(zigm22(mlon0:mlon1,mlat0:mlat1) ,'EDYN_ZIGM22' ,1,mlon0,mlon1,mlat0,mlat1)
975 0 : call savefld_waccm(zigmc (mlon0:mlon1,mlat0:mlat1) ,'EDYN_ZIGMC' ,1,mlon0,mlon1,mlat0,mlat1)
976 0 : call savefld_waccm(zigm2 (mlon0:mlon1,mlat0:mlat1) ,'EDYN_ZIGM2' ,1,mlon0,mlon1,mlat0,mlat1)
977 0 : call savefld_waccm(rim1 (mlon0:mlon1,mlat0:mlat1) ,'EDYN_RIM1' ,1,mlon0,mlon1,mlat0,mlat1)
978 0 : call savefld_waccm(rim2 (mlon0:mlon1,mlat0:mlat1) ,'EDYN_RIM2' ,1,mlon0,mlon1,mlat0,mlat1)
979 : endif
980 :
981 7296 : end subroutine complete_integrals
982 : !-----------------------------------------------------------------------
983 7296 : subroutine rhspde()
984 : use edyn_params, only: pi_dyn, r0
985 : use edyn_maggrid, only: dlatm, dlonm, rcos0s, dt1dts
986 : !
987 : ! Calculate right-hand side from rim1,2 on mag subdomains.
988 : ! Use global longitude arrays for poles and equator obtained
989 : ! by sub complete_integrals.
990 : !
991 : ! Local:
992 : integer :: j,i
993 14592 : real(r8), dimension(nmlat) :: tint1
994 : real(r8) :: &
995 14592 : rim2_npm1(nmlonp1), & ! global rim2 at nmlat-1
996 14592 : rim2_eqp1(nmlonp1), & ! global rim2 at meq+1
997 14592 : rim2_meq (nmlonp1), & ! global rim2 at mag eq
998 14592 : rim2_tmp (nmlonp1), & ! temp array
999 14592 : rim1_meq (nmlonp1), & ! global rim1 at mag eq
1000 14592 : zigm2_meq(nmlonp1), & ! needed for rim1_meq
1001 14592 : zigmc_meq(nmlonp1), & ! needed for rim1_meq
1002 14592 : zigm22_meq(nmlonp1) ! needed for rim1_meq
1003 : real(r8) :: r8_nmlon
1004 7296 : r8_nmlon = real(nmlon, r8)
1005 :
1006 715008 : do j=1,nmlat
1007 715008 : tint1(j) = cos(-pi_dyn/2._r8+(j-1)*dlatm)
1008 : enddo
1009 : !
1010 : ! Init rhs subdomains:
1011 365199 : rhs(:,:) = finit
1012 : !
1013 : ! Will need rim2 at npole-1 and mag equator:
1014 : ! rim2_npm1: global rim2 at nmlat-1:
1015 598272 : rim2_npm1(:) = fpole_jpm2(:,3,6)+fpole_jpm2(:,1,6)
1016 : !
1017 : ! rim2_meq: global rim2 at mag equator:
1018 598272 : rim2_meq(:) = .060_r8*(feq_jpm1(:,1,6)+feq_jpm1(:,2,6))
1019 598272 : rim2_meq(:) = rim2_meq(:)/rcos0s(nmlath)
1020 598272 : rim2_meq(:) = rim2_meq(:)*2._r8 ! fold eq on itself
1021 : !
1022 : ! Perform differentiation of rim(2) w.r.t. lam_0:
1023 : ! +/- (d [ K_(m lam)^D * cos(lam_m)]/ d lam_0 ) /cos ( lam_0) =
1024 : ! + (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) / (RCOS0S*DT0DTS) =
1025 : ! +/- (d [ K_(m lam)^D(0) * cos(lam_0)]/ d lam_0 ) /cos ( lam_0) =
1026 : !
1027 : ! Lat scan to define rhs subdomains:
1028 29412 : do j=mlat0,mlat1
1029 : !
1030 : ! North Pole (redundant in longitude):
1031 29412 : if (j == nmlat) then ! north pole
1032 1767 : do i=mlon0,mlon1
1033 1539 : rhs(i,j) = -2._r8/r8_nmlon*dot_product(unitvm,rim2_npm1(1:nmlon))/ &
1034 126426 : tint1(nmlat-1)
1035 : enddo
1036 : !
1037 : ! Include the boundary condition at the equator.
1038 : ! rhs(equator)/R = d (K_mphi^DT(0) - sig_philam/sig_lamlam*K_mlam^DT(0)) / d phi_m
1039 : ! + d (cos lam_0 * K_mlam^DT(0))/ d lam_0
1040 : ! from Cicely's notes:
1041 : ! I_1 = 0.5*(K_(m phi)^DT(0) - Sig_(phi lam)/Sig_(lam lam)*K_(ml am)^DT(0))
1042 : ! I_2 = K_(m lam)^DT(0)
1043 : ! differentiate
1044 : ! rhs = (I_1(i+1/2,j)-I_1(i-1/2,j))/dlonm +
1045 : ! (2*cos(lam_0)_(j+1/2)*I_2(i,j+1/2))/dlat_0
1046 : ! (first calc global mag equator as in complete_integrals)
1047 : !
1048 21888 : elseif (j == nmlath) then ! mag equator
1049 0 : rim2_eqp1(:) = feq_jpm1(:,2,6)/rcos0s(j-1)+ &
1050 18696 : feq_jpm1(:,1,6)/rcos0s(j+1)
1051 18696 : zigm22_meq(:) = .125_r8*(feq_jpm1(:,1,2)+feq_jpm1(:,2,2))
1052 18696 : zigmc_meq (:) = .125_r8*(feq_jpm1(:,1,3)+feq_jpm1(:,2,3))
1053 18696 : zigm2_meq (:) = .060_r8*(feq_jpm1(:,1,4)+feq_jpm1(:,2,4))
1054 18696 : rim1_meq (:) = .060_r8*(feq_jpm1(:,1,5)+feq_jpm1(:,2,5))
1055 18696 : rim2_tmp (:) = .060_r8*(feq_jpm1(:,1,6)+feq_jpm1(:,2,6))
1056 18696 : do i=1,nmlonp1
1057 18468 : rim1_meq(i) = rim1_meq(i) - (zigm2_meq(i)-zigmc_meq(i))/ &
1058 37164 : zigm22_meq(i)*rim2_tmp(i)
1059 : enddo
1060 18696 : rim1_meq(:) = rim1_meq(:)/dt1dts(j)
1061 18696 : rim1_meq(:) = rim1_meq(:)*2._r8 ! fold eq on itself
1062 :
1063 1767 : do i=mlon0,mlon1
1064 1767 : if (i==1) then ! western most lon
1065 19 : rhs(i,j) = 0.5_r8/dlonm*(rim1(i+1,j)-rim1_meq(nmlon))
1066 19 : rhs(i,j) = rhs(i,j)+1._r8/dlatm*(tint1(j)*rim2(i,j)+ &
1067 38 : tint1(j+1)*rim2_eqp1(i))
1068 1520 : elseif (i==nmlonp1) then ! eastern most lon
1069 19 : rhs(i,j) = 0.5_r8/dlonm*(rim1_meq(1)-rim1(i-1,j))
1070 19 : rhs(i,j) = rhs(i,j)+1._r8/dlatm*(tint1(j)*rim2(i,j)+ &
1071 38 : tint1(j+1)*rim2_eqp1(i))
1072 : else ! body of i subdomain
1073 : ! Note that rim1 halos were set before calling this subroutine.
1074 1501 : rhs(i,j) = 0.5_r8/dlonm*(rim1(i+1,j)-rim1(i-1,j))
1075 1501 : rhs(i,j) = rhs(i,j)+1._r8/dlatm*(tint1(j)*rim2(i,j)+ &
1076 3002 : tint1(j+1)*rim2_eqp1(i))
1077 : endif
1078 : enddo ! i=mlon0,mlon1
1079 : !
1080 : ! North hemisphere (not npole and not equator):
1081 : ! (allow south hemisphere to remain 0)
1082 : ! (use rim1 instead of tint33)
1083 : !
1084 21660 : elseif (j > nmlath) then ! north hem only (excluding npole)
1085 83049 : do i=mlon0,mlon1
1086 72333 : rhs(i,j) = 1._r8/(dlonm*tint1(j))*0.5_r8*(rim1(i+1,j)-rim1(i-1,j))
1087 83049 : if (j == nmlath+1) then
1088 0 : rhs(i,j) = rhs(i,j)+1._r8/(dlatm*tint1(j))* &
1089 1539 : 0.5_r8*(rim2(i,j+1)*tint1(j+1)-rim2_meq(i)*tint1(j-1))
1090 : else
1091 0 : rhs(i,j) = rhs(i,j)+1._r8/(dlatm*tint1(j))* &
1092 70794 : 0.5_r8*(rim2(i,j+1)*tint1(j+1)-rim2(i,j-1)*tint1(j-1))
1093 : endif
1094 : enddo
1095 : endif ! at poles or equator
1096 : enddo ! j=mlat0,mlat1
1097 : !
1098 : ! scale (multiply by earth radius in meter = R0*1.E-2)
1099 : ![( d K_(m phi)^D / d phi /(cos(theta_m)?) +
1100 : ! (d [ K_(m lam)^D * cos(lam_m)]/ d lam_m ) /cos ( lam_m) ] * R / (RCOS0S*DT0DTS)
1101 : ! ~ J_(Mr)*r^2*cos(theta_m)/cos(theta_0)/DT0DTS
1102 : ! theta_m = theta_s
1103 : !
1104 29412 : do j=mlat0,mlat1
1105 178695 : do i=mlon0,mlon1
1106 171399 : rhs(i,j) = rhs(i,j)*r0*1.e-2_r8
1107 : enddo
1108 : enddo
1109 :
1110 7296 : end subroutine rhspde
1111 : !-----------------------------------------------------------------------
1112 7296 : subroutine gather_edyn()
1113 : !
1114 : ! Gather needed global arrays to root task, so it can finish non-parallel
1115 : ! part of dynamo (beginning after sub rhspde) as in original code
1116 : !
1117 : use edyn_mpi, only: mp_gather_edyn
1118 : use edyn_solve, only: & ! (nmlonp1,nmlat)
1119 : zigm11_glb, &
1120 : zigm22_glb, &
1121 : zigmc_glb, &
1122 : zigm2_glb, &
1123 : rhs_glb
1124 : use edyn_solve, only: rim_glb ! pde solver output (nmlonp1,nmlat,2)
1125 : !
1126 : ! Local:
1127 : ! 7 fields to gather: zigm11,zigm22,zigmc,zigm2,rim1,rim2,rhs
1128 : !
1129 : integer, parameter :: nf = 7
1130 14592 : real(r8) :: fmsub(mlon0:mlon1,mlat0:mlat1,nf)
1131 14592 : real(r8) :: fmglb(nmlonp1,nmlat,nf)
1132 14592 : real(r8) :: rhs_nhem(nmlonp1,nmlat)
1133 : integer :: i,j,jj
1134 : !
1135 : ! These calls exclude halo points in zigm11, etc.
1136 : !
1137 185991 : fmsub(:,:,1) = zigm11(mlon0:mlon1,mlat0:mlat1)
1138 178695 : fmsub(:,:,2) = zigm22(mlon0:mlon1,mlat0:mlat1)
1139 178695 : fmsub(:,:,3) = zigmc (mlon0:mlon1,mlat0:mlat1)
1140 178695 : fmsub(:,:,4) = zigm2 (mlon0:mlon1,mlat0:mlat1)
1141 178695 : fmsub(:,:,5) = rim1 (mlon0:mlon1,mlat0:mlat1)
1142 178695 : fmsub(:,:,6) = rim2 (mlon0:mlon1,mlat0:mlat1)
1143 178695 : fmsub(:,:,7) = rhs (mlon0:mlon1,mlat0:mlat1)
1144 :
1145 : call mp_gather_edyn(fmsub,mlon0,mlon1,mlat0,mlat1, &
1146 7296 : fmglb,nmlonp1,nmlat,nf)
1147 : !
1148 : ! Now root task can take over, and work with global arrays:
1149 : !
1150 7296 : if (mytid==0) then
1151 151145 : zigm11_glb(:,:) = fmglb(:,:,1)
1152 151145 : zigm22_glb(:,:) = fmglb(:,:,2)
1153 151145 : zigmc_glb(:,:) = fmglb(:,:,3)
1154 151145 : zigm2_glb(:,:) = fmglb(:,:,4)
1155 151145 : rim_glb(:,:,1) = fmglb(:,:,5)
1156 151145 : rim_glb(:,:,2) = fmglb(:,:,6)
1157 151145 : rhs_nhem(:,:) = fmglb(:,:,7)
1158 : !
1159 : ! Transfer local global rhs_nhem (from sub rhspde) to rhs_glb,
1160 : ! so the latter has data in south hemisphere.
1161 : !
1162 151145 : rhs_glb= 0._r8 ! init
1163 1862 : do j=1,nmlat
1164 : !
1165 : ! Transfer north pole to equator:
1166 1862 : if (j == nmlat) then
1167 1558 : do i=1,nmlonp1
1168 1558 : rhs_glb(i,nmlath) = rhs_nhem(i,j)
1169 : enddo
1170 : ! Transfer equator to south pole:
1171 1824 : elseif (j == nmlath) then
1172 1558 : do i=1,nmlonp1
1173 1558 : rhs_glb(i,1) = rhs_nhem(i,j)
1174 : enddo
1175 : ! Transfer north hem to south hem:
1176 1805 : elseif (j > nmlath) then ! 50 -> 96
1177 893 : jj = j-nmlath+1 ! 2 -> 48
1178 73226 : do i=1,nmlonp1
1179 73226 : rhs_glb(i,jj) = rhs_nhem(i,j)
1180 : enddo
1181 : endif
1182 : enddo ! j=mlat0,mlat1
1183 :
1184 : endif ! mytid==0
1185 7296 : end subroutine gather_edyn
1186 : !-----------------------------------------------------------------------
1187 7296 : subroutine highlat_poten()
1188 : use edyn_solve, only: &
1189 : phihm, & ! high-latitude potential (nmlonp1,nmlat)
1190 : pfrac ! NH fraction of potential (nmlonp1,nmlat0)
1191 : !
1192 : ! Global PDE solution rim_glb(:,:,1) has been scattered to mag subdomains
1193 : ! in rim1, and halos set (this overwrites previous rim1 from fieldline
1194 : ! integrations). Now add high latitude potential from empirical model
1195 : ! (heelis or weimer), defining phim2d on mag subdomains. After this pthreed
1196 : ! will expand phim2d to phim3d.
1197 : !
1198 : ! Input: rim1(mag subdomains) ! Solution from mudpack solver (nhem only)
1199 : ! pfrac(nmlonp1,nmlat0) ! NH fraction of potential
1200 : ! phihm(nmlonp1,nmlat) ! potential in magnetic
1201 : ! Output: phim2d(mag subdomains) ! solution with phihm in both nhem and shem
1202 : ! Both rim1 and phim2d are dimensioned (mlon00:mlon11,mlat00:mlat11)
1203 : !
1204 : ! Both phihm and pfrac have been set by either heelis or weimer.
1205 : ! phihm is on 2d global mag grid, pfrac is in north hemisphere only
1206 : !
1207 : ! Local:
1208 : logical, parameter :: mod_heelis = .false. ! true == modified
1209 : integer :: i,j,jn,js
1210 : real(r8) :: fac
1211 : !
1212 : ! Add empirical model potential at high latitude:
1213 : !
1214 7296 : fac = 1.0_r8
1215 : if (mod_heelis) fac = 0._r8 ! modified heelis
1216 29412 : do j=mlat0,mlat1
1217 22116 : if (j > nmlath) cycle ! south only (including equator)
1218 11172 : jn = nmlat-j+1
1219 11172 : js = nmlath-j+1
1220 93879 : do i=mlon0,mlon1
1221 0 : phim2d(i,j) = rim1(i,j)+fac*(1._r8-pfrac(i,js))*(phihm(i,j)- &
1222 97527 : phihm(i,jn))
1223 : enddo
1224 : enddo
1225 :
1226 29412 : do j=mlat0,mlat1
1227 22116 : if (j <= nmlath) cycle ! north only (excluding equator)
1228 92112 : do i=mlon0,mlon1
1229 95988 : phim2d(i,j) = rim1(i,j)
1230 : enddo
1231 : enddo
1232 :
1233 29412 : do j=mlat0,mlat1
1234 22116 : call outfld('PHIHM',phihm(mlon0:omlon1,j),omlon1-mlon0+1,j)
1235 29412 : call outfld('PHIM2D',phim2d(mlon0:omlon1,j),omlon1-mlon0+1,j)
1236 : enddo
1237 :
1238 7296 : end subroutine highlat_poten
1239 : !-----------------------------------------------------------------------
1240 7296 : subroutine pthreed()
1241 : !
1242 : ! phim2d is now 2d electric potential solution on mag subdomains,
1243 : ! with high-latitude potential added from empirical model (see subs
1244 : ! heelis and highlat_poten), and mag halos set. Now expand phim2d in
1245 : ! vertical, defining phim3d. Also calculate electric field ed13d, ed23d
1246 : ! for later current calculations, and ephi3d, elam3d and emz3d for conversion
1247 : ! to geographic grid (sub pefield), and subsequent calculation of ion drifts
1248 : ! by sub ionvel (not in edynamo).
1249 : !
1250 : use edyn_params, only: Rearth, pi_dyn, r0, kbotdyn
1251 : use edyn_maggrid, only: ylatm, dlatm, dlonm, rcos0s, dt1dts, dt0dts, table
1252 : use edyn_mpi, only: &
1253 : mp_mag_halos, &
1254 : mp_magpole_2d, &
1255 : mp_mageq_jpm3, &
1256 : mp_mag_jslot, &
1257 : mp_magpoles, &
1258 : mp_mag_periodic_f2d, &
1259 : ixfind
1260 : !
1261 : ! Local:
1262 : real(r8), parameter :: eps = 1.e-10_r8
1263 : integer :: mxneed
1264 : integer :: i,j,k,n,mlon00,mlon11,mlat00,mlat11
1265 : real(r8) :: csth0, cosltm, sym, pi, phims, phimn, rind
1266 14592 : real(r8), dimension(nmlonp1) :: thetam,pslot,qslot
1267 14592 : integer, dimension(nmlonp1) :: islot,jslot,ip1f,ip2f,ip3f
1268 :
1269 : ! real(r8), dimension(mlon0-1:mlon1+1,mlat0-1:mlat1+1) :: ed1,ed2
1270 :
1271 14592 : real(r8), dimension(mlon0-1:mlon1+1,mlat0-1:mlat1+1) :: ephi,elam
1272 14592 : real(r8) :: fpole2d_jpm2(nmlonp1,4,4) ! global lons at S pole+1,2 and N pole-1,2
1273 14592 : real(r8) :: fpoles(nmlonp1,2,1) ! global lons at poles (1 field only)
1274 14592 : real(r8) :: fmsub(mlon0:mlon1,mlat0:mlat1,4)
1275 14592 : real(r8) :: fmsub1(mlon0-1:mlon1+1,mlat0-1:mlat1+1,5)
1276 14592 : real(r8) :: feq_jpm3(nmlonp1,-3:3,1) ! global lons at equator +/- 3
1277 14592 : integer :: jneed(nmlat+2) ! lats needed from other tasks for interp
1278 : integer :: njneed,icount
1279 : real(r8), dimension(mlon0-1:mlon1+1,nmlat+2) :: &
1280 14592 : phineed, & ! phim2d at needed latitudes
1281 14592 : ed1need, & ! ed1 at needed latitudes
1282 14592 : ed2need, & ! ed2 at needed latitudes
1283 14592 : ephineed, & ! ephi at needed latitudes
1284 14592 : elamneed ! elam at needed latitudes
1285 14592 : real(r8), dimension(mlon0-1:mlon1+1,nmlat+2,5) :: fmneed
1286 : real(r8) :: phi0j0,phi1j0,phi0j1,phi1j1
1287 : real(r8) :: ed1i0j0,ed1i1j0,ed1i0j1,ed1i1j1
1288 : real(r8) :: ed2i0j0,ed2i1j0,ed2i0j1,ed2i1j1
1289 : real(r8) :: ephi0j0,ephi1j0,ephi0j1,ephi1j1
1290 : real(r8) :: elam0j0,elam1j0,elam0j1,elam1j1
1291 : real(r8) :: fac_elam
1292 : !
1293 7296 : mxneed=nmlat+2
1294 7296 : pi = pi_dyn
1295 7296 : mlon00=mlon0-1 ; mlon11=mlon1+1
1296 7296 : mlat00=mlat0-1 ; mlat11=mlat1+1
1297 : !
1298 : ! Calculate ed1,ed2 components of electric field:
1299 : ! phim2d has halos set, so when mlon0==1, i-1 should wrap
1300 : ! to value at i==nmlon, and when mlon1==nmlonp1, i+1 should
1301 : ! wrap to value at i==2.
1302 : !
1303 365199 : ed1 = 0._r8
1304 365199 : ephi = 0._r8
1305 29412 : do j=mlat0,mlat1
1306 22116 : if (j==1.or.j==nmlat) cycle
1307 21660 : csth0 = cos(-pi/2._r8+(j-1)*dlatm)
1308 175161 : do i=mlon0,mlon1
1309 146205 : ed1(i,j) = -(phim2d(i+1,j)-phim2d(i-1,j))/(2._r8*dlonm*csth0)* &
1310 292410 : rcos0s(j)/(r0*1.e-2_r8)
1311 168321 : ephi(i,j) = ed1(i,j)*(r0*1.e-2_r8)
1312 : enddo ! i=mlon0,mlon1
1313 : enddo ! j=mlat0,mlat1
1314 : !
1315 : ! Southern hemisphere (excluding equator):
1316 365199 : ed2 = 0._r8
1317 365199 : elam = 0._r8
1318 29412 : do j=mlat0,mlat1
1319 22116 : if (j >= nmlath) cycle
1320 10944 : if (j==1.or.j==nmlat) cycle ! skip poles
1321 90345 : do i=mlon0,mlon1
1322 72333 : ed2(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2._r8*dlatm)*dt1dts(j)/ &
1323 144666 : (r0*1.e-2_r8)
1324 94449 : elam(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2._r8*dlatm)*dt0dts(j)
1325 : enddo ! i=mlon0,mlon1
1326 : enddo ! j=mlat0,mlat1
1327 : !
1328 : ! Northern hemisphere (excluding equator):
1329 29412 : do j=mlat0,mlat1
1330 22116 : if (j <= nmlath) cycle
1331 10944 : if (j==1.or.j==nmlat) cycle ! skip poles
1332 90345 : do i=mlon0,mlon1
1333 72333 : ed2(i,j) = (phim2d(i,j+1)-phim2d(i,j-1))/(2._r8*dlatm)*dt1dts(j)/ &
1334 144666 : (r0*1.e-2_r8)
1335 94449 : elam(i,j)= -(phim2d(i,j+1)-phim2d(i,j-1))/(2._r8*dlatm)*dt0dts(j)
1336 : enddo ! i=mlon0,mlon1
1337 : enddo ! j=mlat0,mlat1
1338 :
1339 : ! Need ed1,2 at global longitudes at j==2 and j==nmlat-1:
1340 : ! mp_magpole_2d: Return fpole_jpm2(nglblon,1->4,nf) as:
1341 : ! 1: j = jspole+1 (spole+1)
1342 : ! 2: j = jspole+2 (spole+2) (unused here)
1343 : ! 3: j = jnpole-1 (npole-1)
1344 : ! 4: j = jnpole-2 (npole-2) (unused here)
1345 : !
1346 178695 : fmsub(:,:,1) = ed1(mlon0:mlon1,mlat0:mlat1)
1347 178695 : fmsub(:,:,2) = ed2(mlon0:mlon1,mlat0:mlat1)
1348 178695 : fmsub(:,:,3) = ephi(mlon0:mlon1,mlat0:mlat1)
1349 178695 : fmsub(:,:,4) = elam(mlon0:mlon1,mlat0:mlat1)
1350 : call mp_magpole_2d(fmsub,mlon0,mlon1,mlat0,mlat1,nmlonp1, &
1351 7296 : 1,nmlat,fpole2d_jpm2,4)
1352 : !
1353 : ! Poles: average over four surrounding points
1354 598272 : do i = 1,nmlonp1
1355 590976 : ip1f(i) = i + nmlon/4
1356 590976 : if (ip1f(i) > nmlonp1) ip1f(i) = ip1f(i) - nmlon
1357 590976 : ip2f(i) = i + nmlon/2
1358 590976 : if (ip2f(i) > nmlonp1) ip2f(i) = ip2f(i) - nmlon
1359 590976 : ip3f(i) = i + 3*nmlon/4
1360 598272 : if (ip3f(i) > nmlonp1) ip3f(i) = ip3f(i) - nmlon
1361 : enddo ! i=1,nmlonp1
1362 : !
1363 : ! S pole:
1364 29412 : do j=mlat0,mlat1
1365 29412 : if (j==1) then
1366 1767 : do i=mlon0,mlon1
1367 3078 : ed1(i,j)=.25_r8*(fpole2d_jpm2(i,1,1)-fpole2d_jpm2(ip2f(i),1,1)+ &
1368 4617 : fpole2d_jpm2(ip1f(i),1,2)-fpole2d_jpm2(ip3f(i),1,2))
1369 1539 : ed2(i,j)=.25_r8*(fpole2d_jpm2(i,1,2)-fpole2d_jpm2(ip2f(i),1,2)- &
1370 3078 : fpole2d_jpm2(ip1f(i),1,1)+fpole2d_jpm2(ip3f(i),1,1))
1371 :
1372 3078 : ephi(i,j)=.25_r8*(fpole2d_jpm2(i,1,3)-fpole2d_jpm2(ip2f(i),1,3)+ &
1373 3078 : fpole2d_jpm2(ip1f(i),1,4)-fpole2d_jpm2(ip3f(i),1,4))
1374 1539 : elam(i,j)=.25_r8*(fpole2d_jpm2(i,1,4)-fpole2d_jpm2(ip2f(i),1,4)- &
1375 3306 : fpole2d_jpm2(ip1f(i),1,3)+fpole2d_jpm2(ip3f(i),1,3))
1376 : enddo ! i=mlon0,mlon1
1377 : ! N pole:
1378 21888 : elseif (j==nmlat) then
1379 1767 : do i=mlon0,mlon1
1380 3078 : ed1(i,j)=.25_r8*(fpole2d_jpm2(i,3,1)-fpole2d_jpm2(ip2f(i),3,1)+ &
1381 4617 : fpole2d_jpm2(ip1f(i),3,2)-fpole2d_jpm2(ip3f(i),3,2))
1382 1539 : ed2(i,j)=.25_r8*(fpole2d_jpm2(i,3,2)-fpole2d_jpm2(ip2f(i),3,2)- &
1383 3078 : fpole2d_jpm2(ip1f(i),3,1)+fpole2d_jpm2(ip3f(i),3,1))
1384 :
1385 3078 : ephi(i,j)=.25_r8*(fpole2d_jpm2(i,3,3)-fpole2d_jpm2(ip2f(i),3,3)+ &
1386 3078 : fpole2d_jpm2(ip1f(i),3,4)-fpole2d_jpm2(ip3f(i),3,4))
1387 1539 : elam(i,j)=.25_r8*(fpole2d_jpm2(i,3,4)-fpole2d_jpm2(ip2f(i),3,4)- &
1388 3306 : fpole2d_jpm2(ip1f(i),3,3)+fpole2d_jpm2(ip3f(i),3,3))
1389 : enddo ! i=mlon0,mlon1
1390 : endif ! S or N pole
1391 : enddo ! j=mlat0,mlat1
1392 : !
1393 : ! Equator: derivative of quadratic polynomial (3 given points)
1394 : ! For equator and equator +/- 1 of ed2, we need equator and
1395 : ! equator +/- 3 of phim2d (note feq_jpm3(nmlonp1,-3:3,1)):
1396 : !
1397 178695 : fmsub(:,:,1) = phim2d(mlon0:mlon1,mlat0:mlat1)
1398 7296 : call mp_mageq_jpm3(fmsub(:,:,1),mlon0,mlon1,mlat0,mlat1,nmlonp1,feq_jpm3,1)
1399 29412 : do j=mlat0,mlat1
1400 29412 : if (j==nmlath-1) then ! equator-1
1401 1767 : do i=mlon0,mlon1
1402 1539 : ed2(i,j) = (4._r8*feq_jpm3(i,-2,1)-feq_jpm3(i,-3,1)- &
1403 3306 : 3._r8*feq_jpm3(i,-1,1))/(2._r8*dlatm)/(r0*1.e-2_r8)
1404 : enddo
1405 21888 : elseif (j==nmlath) then ! equator
1406 1767 : do i=mlon0,mlon1
1407 1539 : ed2(i,j) = (4._r8*feq_jpm3(i,1,1)-feq_jpm3(i,2,1)- &
1408 1539 : 3._r8*feq_jpm3(i,0,1))/(2._r8*dlatm)/(r0*1.e-2_r8)
1409 1539 : elam(i,j) = (4._r8*feq_jpm3(i,1,1)-feq_jpm3(i,2,1)- &
1410 1767 : 3._r8*feq_jpm3(i,0,1))/(2._r8*dlatm)
1411 : enddo
1412 21660 : elseif (j==nmlath+1) then ! equator+1
1413 1767 : do i=mlon0,mlon1
1414 1539 : ed2(i,j) = (4._r8*feq_jpm3(i,2,1)-feq_jpm3(i,3,1)- &
1415 3306 : 3._r8*feq_jpm3(i,1,1))/(2._r8*dlatm)/(r0*1.e-2_r8)
1416 : enddo
1417 : endif ! equator +/- 1
1418 : enddo ! j=mlat0,mlat1
1419 : !
1420 : ! Set halos for 3d calculations:
1421 365199 : fmsub1(:,:,1) = ed1
1422 365199 : fmsub1(:,:,2) = ed2
1423 365199 : fmsub1(:,:,3) = ephi
1424 365199 : fmsub1(:,:,4) = elam
1425 7296 : call mp_mag_halos(fmsub1,mlon0,mlon1,mlat0,mlat1,4)
1426 372495 : ed1 = fmsub1(:,:,1)
1427 372495 : ed2 = fmsub1(:,:,2)
1428 365199 : ephi = fmsub1(:,:,3)
1429 365199 : elam = fmsub1(:,:,4)
1430 :
1431 29412 : do j=mlat0,mlat1
1432 22116 : call outfld('ED1',ed1(mlon0:omlon1,j),omlon1-mlon0+1,j)
1433 29412 : call outfld('ED2',ed2(mlon0:omlon1,j),omlon1-mlon0+1,j)
1434 : enddo
1435 : !
1436 : ! Determine latitudes needed for interpolation that fall
1437 : ! outside a task's latitudinal subdomain:
1438 : !
1439 7296 : if (debug) write(iulog,*) "pthreed: kbotdyn ", kbotdyn
1440 :
1441 7296 : njneed = 0 ! number of unique latitudes needed
1442 729600 : jneed(:) = -1 ! j-indices of needed latitudes
1443 510720 : do k=kbotdyn,nmlev
1444 2036724 : do j=mlat0,mlat1
1445 1526004 : if (j==1.or.j==nmlat) cycle ! exclude poles
1446 1494540 : sym = 1._r8
1447 1494540 : if (j < nmlath) sym = -1._r8
1448 1494540 : cosltm = cos(ylatm(j))
1449 12086109 : do i=mlon0,mlon1
1450 10088145 : if (i==nmlonp1) cycle
1451 :
1452 9963600 : thetam(i)=(Rearth+zpotm3d(i,j,kbotdyn))/(Rearth+zpotm3d(i,j,k))
1453 9963600 : thetam(i) = acos(sqrt(thetam(i))*cosltm*(1._r8-eps))
1454 :
1455 9963600 : pslot(i) = thetam(i)*180._r8/pi+1._r8
1456 9963600 : islot(i) = pslot(i)
1457 9963600 : rind = real(islot(i), kind=r8)
1458 9963600 : pslot(i) = pslot(i)-rind
1459 :
1460 9963600 : thetam(i) = ((1._r8-pslot(i))*table(islot(i),2)+pslot(i)* &
1461 19927200 : table(islot(i)+1,2))*sym ! thetam negative for south hem
1462 :
1463 9963600 : islot(i) = i
1464 9963600 : pslot(i) = 0._r8
1465 9963600 : qslot(i) = (thetam(i)+pi/2._r8)/dlatm+1._r8
1466 9963600 : jslot(i) = qslot(i)
1467 9963600 : rind = real(jslot(i), kind=r8)
1468 9963600 : qslot(i) = qslot(i)-rind
1469 :
1470 : ! Save j index if outside subdomain w/ halos:
1471 753924862 : if ((jslot(i) < mlat00 .or. jslot(i) > mlat11).and. &
1472 : .not.any(jslot(i)==jneed)) then
1473 14146 : njneed = njneed+1
1474 14146 : if (njneed > mxneed) call endrun('njneed')
1475 14146 : jneed(njneed) = jslot(i)
1476 : endif ! jslot is outside subdomain
1477 : !
1478 : ! Save j+1 index if outside subdomain:
1479 754239999 : if ((jslot(i)+1 < mlat00 .or. jslot(i)+1 > mlat11).and. &
1480 1526004 : .not.any(jslot(i)+1==jneed)) then
1481 17052 : njneed = njneed+1
1482 17052 : if (njneed > mxneed) call endrun('njneed')
1483 17052 : jneed(njneed) = jslot(i)+1
1484 : endif ! jslot(i)+1 is outside subdomain
1485 : enddo ! i=mlon0,mlon1
1486 : enddo ! j=mlat0,mlat1
1487 : enddo ! k=kbotdyn,nmlev
1488 : !
1489 : ! Get phim2 at needed latitudes (note inclusion of phim2d halos).
1490 : ! real,intent(in) :: fin(mlon00:mlon11,mlat00:mlat11,nf) ! data at current subdomain
1491 : ! real,intent(out) :: fout(mlon00:mlon11,mxneed,nf) ! returned data at needed lats
1492 : !
1493 365199 : fmsub1(:,:,1) = phim2d
1494 365199 : fmsub1(:,:,2) = ed1
1495 365199 : fmsub1(:,:,3) = ed2
1496 365199 : fmsub1(:,:,4) = ephi
1497 365199 : fmsub1(:,:,5) = elam
1498 7296 : call mp_mag_jslot(fmsub1,mlon00,mlon11,mlat00,mlat11,fmneed,jneed,mxneed,5)
1499 7049760 : phineed = fmneed(:,:,1)
1500 7049760 : ed1need = fmneed(:,:,2)
1501 7049760 : ed2need = fmneed(:,:,3)
1502 7049760 : ephineed= fmneed(:,:,4)
1503 7049760 : elamneed= fmneed(:,:,5)
1504 :
1505 23237646 : ephi3d = 0._r8
1506 23237646 : elam3d = 0._r8
1507 23237646 : emz3d = 0._r8
1508 510720 : do k=kbotdyn,nmlev
1509 2036724 : do j=mlat0,mlat1
1510 1526004 : if (j==1.or.j==nmlat) cycle ! exclude poles
1511 1494540 : sym = 1._r8
1512 1494540 : if (j < nmlath) sym = -1._r8
1513 1494540 : cosltm = cos(ylatm(j))
1514 12086109 : do i=mlon0,mlon1
1515 10088145 : if (i==nmlonp1) cycle
1516 :
1517 9963600 : thetam(i)=(Rearth+zpotm3d(i,j,kbotdyn))/(Rearth+zpotm3d(i,j,k))
1518 9963600 : thetam(i) = acos(sqrt(thetam(i))*cosltm*(1._r8-eps))
1519 9963600 : fac_elam = tan(ylatm(j))/tan(thetam(i)*sym) ! tan(lambda_q)/tan(lambda_m)
1520 :
1521 9963600 : pslot(i) = thetam(i)*180._r8/pi+1._r8
1522 9963600 : islot(i) = pslot(i)
1523 9963600 : rind = real(islot(i), kind=r8)
1524 9963600 : pslot(i) = pslot(i)-rind
1525 :
1526 9963600 : thetam(i) = ((1._r8-pslot(i))*table(islot(i),2)+pslot(i)* &
1527 19927200 : table(islot(i)+1,2))*sym ! thetam negative for south hem
1528 :
1529 9963600 : islot(i) = i
1530 9963600 : pslot(i) = 0._r8
1531 9963600 : qslot(i) = (thetam(i)+pi/2._r8)/dlatm+1._r8
1532 9963600 : jslot(i) = qslot(i)
1533 9963600 : rind = real(jslot(i), kind=r8)
1534 9963600 : qslot(i) = qslot(i)-rind
1535 : !
1536 : ! Check for jslot in subdomain:
1537 9963600 : if (jslot(i) >= mlat00.and.jslot(i) <= mlat11) then ! within subdomain
1538 7341826 : phi0j0 = phim2d(islot(i) ,jslot(i))
1539 7341826 : phi1j0 = phim2d(islot(i)+1,jslot(i))
1540 7341826 : ed1i0j0 = ed1(islot(i) ,jslot(i))
1541 7341826 : ed1i1j0 = ed1(islot(i)+1,jslot(i))
1542 7341826 : ed2i0j0 = ed2(islot(i) ,jslot(i))
1543 7341826 : ed2i1j0 = ed2(islot(i)+1,jslot(i))
1544 7341826 : ephi0j0 = ephi(islot(i) ,jslot(i))
1545 7341826 : ephi1j0 = ephi(islot(i)+1,jslot(i))
1546 7341826 : elam0j0 = elam(islot(i) ,jslot(i))
1547 7341826 : elam1j0 = elam(islot(i)+1,jslot(i))
1548 : else ! jslot outside subdomain
1549 2621774 : n = ixfind(jneed,mxneed,jslot(i),icount)
1550 2621774 : if (n==0) then
1551 0 : write(iulog,"('>>> pthreed: i=',i4,' j=',i4,' k=',i4)") i,j,k
1552 : write(iulog,"(' Could not find jslot ',i4,' in jneed=',/,(10i4))") &
1553 0 : i,j,k,jslot(i),jneed
1554 0 : call endrun('jslot(i) not in jneed')
1555 : endif
1556 2621774 : phi0j0 = phineed(islot(i) ,n)
1557 2621774 : phi1j0 = phineed(islot(i)+1,n)
1558 2621774 : ed1i0j0 = ed1need(islot(i) ,n)
1559 2621774 : ed1i1j0 = ed1need(islot(i)+1,n)
1560 2621774 : ed2i0j0 = ed2need(islot(i) ,n)
1561 2621774 : ed2i1j0 = ed2need(islot(i)+1,n)
1562 2621774 : ephi0j0 =ephineed(islot(i) ,n)
1563 2621774 : ephi1j0 =ephineed(islot(i)+1,n)
1564 2621774 : elam0j0 =elamneed(islot(i) ,n)
1565 2621774 : elam1j0 =elamneed(islot(i)+1,n)
1566 : endif
1567 : !
1568 : ! Check for jslot+1 in subdomain:
1569 9963600 : if (jslot(i)+1 >= mlat00.and.jslot(i)+1 <= mlat11) then ! within subdomain
1570 7342633 : phi0j1 = phim2d(islot(i) ,jslot(i)+1)
1571 7342633 : phi1j1 = phim2d(islot(i)+1,jslot(i)+1)
1572 7342633 : ed1i0j1 = ed1(islot(i) ,jslot(i)+1)
1573 7342633 : ed1i1j1 = ed1(islot(i)+1,jslot(i)+1)
1574 7342633 : ed2i0j1 = ed2(islot(i) ,jslot(i)+1)
1575 7342633 : ed2i1j1 = ed2(islot(i)+1,jslot(i)+1)
1576 7342633 : ephi0j1 = ephi(islot(i) ,jslot(i)+1)
1577 7342633 : ephi1j1 = ephi(islot(i)+1,jslot(i)+1)
1578 7342633 : elam0j1 = elam(islot(i) ,jslot(i)+1)
1579 7342633 : elam1j1 = elam(islot(i)+1,jslot(i)+1)
1580 : else ! jslot+1 outside subdomain
1581 2620967 : n = ixfind(jneed,mxneed,jslot(i)+1,icount)
1582 2620967 : if (n==0) then
1583 0 : write(iulog,"('>>> pthreed: i=',i4,' j=',i4,' k=',i4)") i,j,k
1584 : write(iulog,"(' Could not find jslot+1 ',i4,' in jneed=',/,(10i4))") &
1585 0 : i,j,k,jslot(i)+1,jneed
1586 0 : call endrun('jslot(i)+1 not in jneed')
1587 : endif
1588 2620967 : phi0j1 = phineed(islot(i) ,n)
1589 2620967 : phi1j1 = phineed(islot(i)+1,n)
1590 2620967 : ed1i0j1 = ed1need(islot(i) ,n)
1591 2620967 : ed1i1j1 = ed1need(islot(i)+1,n)
1592 2620967 : ed2i0j1 = ed2need(islot(i) ,n)
1593 2620967 : ed2i1j1 = ed2need(islot(i)+1,n)
1594 2620967 : ephi0j1 =ephineed(islot(i) ,n)
1595 2620967 : ephi1j1 =ephineed(islot(i)+1,n)
1596 2620967 : elam0j1 =elamneed(islot(i) ,n)
1597 2620967 : elam1j1 =elamneed(islot(i)+1,n)
1598 : endif
1599 : !
1600 : ! Do the interpolation:
1601 9963600 : phim3d(i,j,k) = (1._r8-qslot(i))*((1._r8-pslot(i))* &
1602 : phi0j0 + pslot(i) * phi1j0) + qslot(i)*((1._r8-pslot(i))* &
1603 9963600 : phi0j1 + pslot(i) * phi1j1)
1604 :
1605 0 : ed13d(i,j,k) = (1._r8-qslot(i))*((1._r8-pslot(i))* &
1606 : ed1i0j0 + pslot(i) * ed1i1j0) + qslot(i)*((1._r8-pslot(i))* &
1607 9963600 : ed1i0j1 + pslot(i) * ed1i1j1)
1608 0 : ed23d(i,j,k) = (1._r8-qslot(i))*((1._r8-pslot(i))* &
1609 : ed2i0j0 + pslot(i) * ed2i1j0) + qslot(i)*((1._r8-pslot(i))* &
1610 9963600 : ed2i0j1 + pslot(i) * ed2i1j1)
1611 :
1612 0 : ephi3d(i,j,k) = (1._r8-qslot(i))*((1._r8-pslot(i))* &
1613 : ephi0j0 + pslot(i) * ephi1j0) + qslot(i)*((1._r8-pslot(i))* &
1614 9963600 : ephi0j1 + pslot(i) * ephi1j1)
1615 0 : elam3d(i,j,k) = (1._r8-qslot(i))*((1._r8-pslot(i))* &
1616 : elam0j0 + pslot(i) * elam1j0) + qslot(i)*((1._r8-pslot(i))* &
1617 9963600 : elam0j1 + pslot(i) * elam1j1)
1618 11614149 : elam3d(i,j,k) = elam3d(i,j,k)*fac_elam ! add height variation
1619 : enddo ! i=mlon0,mlon1
1620 : enddo ! j=mlat0,mlat1
1621 : enddo ! k=kbotdyn,nmlev
1622 :
1623 : !
1624 : ! Mag poles for phim:
1625 : ! mp_magpoles returns global longitudes at S,N poles in fpoles(nglblon,2,nf)
1626 : !
1627 0 : call mp_magpoles(phim2d(mlon0:mlon1,mlat0:mlat1), &
1628 178695 : mlon0,mlon1,mlat0,mlat1,nmlonp1,1,nmlat,fpoles,1)
1629 :
1630 7296 : rind = real(nmlon, kind=r8)
1631 590976 : phims=dot_product(unitvm,fpoles(1:nmlon,1,1))/rind
1632 590976 : phimn=dot_product(unitvm,fpoles(1:nmlon,2,1))/rind
1633 :
1634 510720 : do k=kbotdyn,nmlev
1635 2036724 : do j=mlat0,mlat1
1636 2029428 : if (j==1) then
1637 121923 : do i=mlon0,mlon1
1638 106191 : phim3d(i,j,k) = phims
1639 106191 : ed13d(i,j,k) = ed1(i,j)
1640 106191 : ed23d(i,j,k) = ed2(i,j)
1641 106191 : ephi3d(i,j,k) = ephi(i,j)
1642 121923 : elam3d(i,j,k) = ed2(i,j)*(r0*1.e-2_r8)
1643 : enddo
1644 1510272 : elseif (j==nmlat) then
1645 121923 : do i=mlon0,mlon1
1646 106191 : phim3d(i,j,k) = phimn
1647 106191 : ed13d(i,j,k) = ed1(i,j)
1648 106191 : ed23d(i,j,k) = ed2(i,j)
1649 106191 : ephi3d(i,j,k) = ephi(i,j)
1650 121923 : elam3d(i,j,k) = -ed2(i,j)*(r0*1.e-2_r8)
1651 : enddo
1652 : endif ! poles
1653 : enddo ! j=mlat0,mlat1
1654 : enddo ! k=kbotdyn,nmlev
1655 : !
1656 : ! Extend kbotdyn downward redundantly:
1657 452352 : do k=1,kbotdyn-1
1658 10900395 : phim3d(:,:,k) = phim3d(:,:,kbotdyn)
1659 10900395 : ephi3d(:,:,k) = ephi3d(:,:,kbotdyn)
1660 10907691 : elam3d(:,:,k) = elam3d(:,:,kbotdyn)
1661 : enddo
1662 : !
1663 : ! Upward electric field:
1664 503424 : do k=kbotdyn,nmlev-1
1665 2007312 : do j=mlat0,mlat1
1666 12151260 : do i=mlon0,mlon1
1667 11655132 : emz3d(i,j,k) = -(phim3d(i,j,k+1)-phim3d(i,j,k-1))
1668 : enddo
1669 : enddo
1670 : enddo
1671 : !
1672 955776 : do k=mlev0,mlev1
1673 955776 : call mp_mag_periodic_f2d(phim3d(:,:,k),mlon0,mlon1,mlat0,mlat1,1)
1674 : enddo
1675 : !
1676 7296 : if (debug_hist) then
1677 0 : do j=mlat0,mlat1
1678 0 : call outfld('EPHI3D',ephi3d(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1679 0 : call outfld('ELAM3D',elam3d(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1680 0 : call outfld('EMZ3D', emz3d(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1681 : enddo
1682 : endif
1683 :
1684 29412 : do j=mlat0,mlat1
1685 22064396 : call outfld('PHIM3D',phim3d(mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1686 22064396 : call outfld('ED13D' ,ed13d (mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1687 22071692 : call outfld('ED23D' ,ed23d (mlon0:omlon1,j,mlev1:mlev0:-1),omlon1-mlon0+1,j)
1688 : enddo
1689 :
1690 7296 : end subroutine pthreed
1691 : !-----------------------------------------------------------------------
1692 7296 : subroutine pefield()
1693 : use edyn_params, only: pi
1694 : use edyn_maggrid, only: dt0dts, dlatm, dlonm, rcos0s
1695 : use edyn_geogrid, only: nlev
1696 : use edyn_mpi, only: mp_magpole_3d, mp_mag_halos, mp_magpoles
1697 : use regridder, only: regrid_mag2geo_3d
1698 : !
1699 : ! Local:
1700 : integer :: i, ii, j, k
1701 : real(r8) :: &
1702 14592 : phi3d(mlon0-1:mlon1+1,mlat0-1:mlat1+1,nmlev), & ! local phi w/ halos
1703 14592 : fpole3d_jpm2(nmlonp1,4,nmlev,1) ! global lons at S pole+1,2 and N pole-1,2
1704 : real(r8) :: csth0
1705 14592 : real(r8) :: fpoles(nmlonp1,2,nmlev) ! global lons at poles
1706 14592 : real(r8), dimension(lon0:lon1,lat0:lat1,nlev) :: exgeo, eygeo, ezgeo
1707 :
1708 : !
1709 : ! Copy phim3d to local phi3d, and set halo points:
1710 29412 : do j = mlat0, mlat1
1711 178695 : do i = mlon0, mlon1
1712 19578189 : phi3d(i,j,:) = phim3d(i,j,:)
1713 : end do
1714 : end do
1715 7296 : call mp_mag_halos(phi3d, mlon0, mlon1, mlat0, mlat1, nmlev)
1716 : !
1717 : ! Return fpole3d_jpm2(nglblon,1->4,nlev,nf) as:
1718 : ! 1: j = jspole+1 (spole+1)
1719 : ! 2: j = jspole+2 (spole+2) not used here
1720 : ! 3: j = jnpole-1 (npole-1)
1721 : ! 4: j = jnpole-2 (npole-2) not used here
1722 : !
1723 0 : call mp_magpole_3d(phim3d(mlon0:mlon1,mlat0:mlat1,:), mlon0, &
1724 7296 : mlon1, mlat0, mlat1, nmlev, nmlonp1, 1, nmlat, fpole3d_jpm2, 1)
1725 : !
1726 : ! Set j=0 and j=nmlat+1 of local phi3d. This overwrites the far
1727 : ! north and south halo points set by mp_mag_halos above.
1728 29412 : do j = mlat0, mlat1
1729 29412 : if (j==1) then
1730 1767 : do i = mlon0, mlon1
1731 1539 : ii = 1 + mod(i-1+nmlon/2,nmlon) ! over the south pole
1732 201837 : phi3d(i,j-1,:) = fpole3d_jpm2(ii,1,:,1)
1733 : end do
1734 21888 : else if (j == nmlat) then
1735 1767 : do i = mlon0, mlon1
1736 1539 : ii = 1 + mod(i-1+nmlon/2,nmlon) ! over the north pole
1737 201837 : phi3d(i,j+1,:) = fpole3d_jpm2(ii,3,:,1)
1738 : end do
1739 : end if ! poles or not
1740 : end do ! j=mlat0,mlat1
1741 : !
1742 : ! Meridional component of electric field:
1743 29412 : do j = mlat0, mlat1
1744 178695 : do i = mlon0, mlon1
1745 447849 : elam3d(i,j,:) = -(phi3d(i,j+1,:)-phi3d(i,j-1,:)) / &
1746 20026038 : (2._r8*dlatm)*dt0dts(j)
1747 : end do
1748 : end do
1749 : !
1750 : ! Zonal component of electric field:
1751 29412 : do j = mlat0, mlat1
1752 22116 : if (j==1 .or. j==nmlat) cycle
1753 21660 : csth0 = cos((-pi / 2._r8) + (real(j-1,kind=r8) * dlatm))
1754 175161 : do i = mlon0, mlon1
1755 438615 : ephi3d(i,j,:) = -(phi3d(i+1,j,:) - phi3d(i-1,j,:)) / &
1756 19613586 : (2._r8 * dlonm * csth0) * rcos0s(j)
1757 : end do
1758 : end do
1759 : !
1760 : ! Polar values for ephi3d (need global lons at poles of elam3d):
1761 7296 : call mp_magpoles(elam3d,mlon0,mlon1,mlat0,mlat1,nmlonp1,1,nmlat,fpoles,nmlev)
1762 29412 : do j = mlat0, mlat1
1763 29412 : if (j == 1) then ! south pole
1764 1767 : do i = mlon0, mlon1
1765 1539 : ii = 1 + mod(i-1+(nmlon/4),nmlon) ! over the south pole
1766 201837 : ephi3d(i,j,:) = fpoles(ii,1,:)
1767 : end do
1768 21888 : else if (j == nmlat) then ! north pole
1769 1767 : do i = mlon0, mlon1
1770 1539 : ii = 1+mod(i-1+((3*nmlon)/4),nmlon) ! over the north pole
1771 201837 : ephi3d(i,j,:) = fpoles(ii,2,:)
1772 : end do
1773 : end if ! poles or not
1774 : end do ! j=mlat0,mlat1
1775 : !
1776 : ! emz = d(phi)/dz
1777 941184 : do k = 2, nmlev-1
1778 3772032 : do j = mlat0, mlat1
1779 22872960 : do i = mlon0, mlon1
1780 21939072 : emz3d(i,j,k) = -(phim3d(i,j,k+1)-phi3d(i,j,k-1))
1781 : end do
1782 : end do
1783 : end do ! k=2,nmlev-1
1784 :
1785 : ! regrid from mag grid to geo grid
1786 7296 : call regrid_mag2geo_3d( ephi3d, exgeo )
1787 7296 : call regrid_mag2geo_3d( elam3d, eygeo )
1788 7296 : call regrid_mag2geo_3d( emz3d, ezgeo )
1789 7296 : call regrid_mag2geo_3d( phim3d, phig3d )
1790 : !
1791 : ! Define ex,ey,ez on geographic subdomains for ionvel:
1792 29184 : do j = lat0, lat1
1793 291840 : do i = lon0, lon1
1794 34407936 : ex(:,i,j) = exgeo(i,j,:)
1795 34407936 : ey(:,i,j) = eygeo(i,j,:)
1796 34407936 : ez(:,i,j) = ezgeo(i,j,:)
1797 34429824 : poten(:,i,j) = phig3d(i,j,:)
1798 : end do
1799 : end do
1800 :
1801 : ! ex,ey,ez(nlev,lon0-2,lon1+2,lat0:lat1)
1802 7296 : if (debug) then
1803 : write(iulog,"(a,2e12.4,' ey=',2e12.4,' ez=',2e12.4)") &
1804 0 : 'pefield after mag2phys: ex=', &
1805 0 : minval(ex(:,lon0:lon1,:)),maxval(ex(:,lon0:lon1,:)), &
1806 0 : minval(ey(:,lon0:lon1,:)),maxval(ey(:,lon0:lon1,:)), &
1807 0 : minval(ez(:,lon0:lon1,:)),maxval(ez(:,lon0:lon1,:))
1808 : end if
1809 :
1810 7296 : call savefld_waccm(poten(1:nlev,lon0:lon1,lat0:lat1),'POTEN', &
1811 14592 : nlev,lon0,lon1,lat0,lat1)
1812 :
1813 7296 : end subroutine pefield
1814 : !-----------------------------------------------------------------------
1815 : !-----------------------------------------------------------------------
1816 :
1817 7296 : subroutine ionvel(z,ui,vi,wi,lon0,lon1, lat0,lat1, lev0,lev1)
1818 : !
1819 : ! Calculate 3d ExB ion drifts from electric field (sub pefield)
1820 : ! on geographic grid.
1821 : !
1822 7296 : use edyn_params, only: Rearth
1823 : use edyn_geogrid, only: nlev
1824 : use getapex, only: rjac ! (nlon+1,jspole:jnpole,2,2)
1825 : use getapex, only: bmod ! magnitude of mag field (nlon+1,jspole:jnpole)
1826 : use getapex, only: xb,yb,zb ! north,east,down mag field (nlon+1,jspole:jnpole)
1827 : !
1828 : ! Args:
1829 : integer,intent(in) :: & ! geographic subdomain
1830 : lon0, lon1, & ! first,last longitude indices of geographic subdomain
1831 : lat0, lat1, & ! first,last latitude indices of geographic subdomain
1832 : lev0, lev1 ! first,last level indices (not distributed)
1833 : real(r8),intent(in), dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: &
1834 : z ! geopotential from input (cm)
1835 : real(r8),intent(out), dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: &
1836 : ui,vi,wi
1837 : !
1838 : ! Local:
1839 : integer :: i,k,j
1840 14592 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: eex,eey,eez
1841 14592 : real(r8), dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: rjac_out
1842 :
1843 : ! mag field diagnostics
1844 291840 : call savefld_waccm(bmod(lon0:lon1,lat0:lat1),'BMOD',1,lon0,lon1,lat0,lat1)
1845 291840 : call savefld_waccm(xb(lon0:lon1,lat0:lat1),'XB',1,lon0,lon1,lat0,lat1)
1846 291840 : call savefld_waccm(yb(lon0:lon1,lat0:lat1),'YB',1,lon0,lon1,lat0,lat1)
1847 291840 : call savefld_waccm(zb(lon0:lon1,lat0:lat1),'ZB',1,lon0,lon1,lat0,lat1)
1848 :
1849 : !
1850 : ! Scan geographic latitude subdomain:
1851 : !
1852 29184 : do j=lat0,lat1
1853 284544 : do i=lon0,lon1
1854 34429824 : do k=lev0,lev1
1855 34145280 : eex(k,i,j) = (rjac(i,j,1,1)*ex(k,i,j)+rjac(i,j,2,1)*ey(k,i,j))/(Rearth+z(k,i,j)) ! V/cm
1856 34407936 : eey(k,i,j) = (rjac(i,j,1,2)*ex(k,i,j)+rjac(i,j,2,2)*ey(k,i,j))/(Rearth+z(k,i,j))
1857 : enddo
1858 : enddo
1859 :
1860 284544 : do i=lon0,lon1
1861 33904512 : do k=lev0+1,lev1-1
1862 33882624 : eez(k,i,j) = ez(k,i,j)/(z(k+1,i,j)-z(k-1,i,j))
1863 : enddo
1864 : enddo
1865 : !
1866 : ! Extrapolate for lower and upper boundaries:
1867 284544 : do i=lon0,lon1
1868 262656 : eez(lev0,i,j) = 2._r8*eez(2,i,j)-eez(3,i,j)
1869 284544 : eez(lev1,i,j) = 2._r8*eez(lev1-1,i,j)-eez(lev1-2,i,j)
1870 : enddo
1871 : !
1872 : ! ion velocities = (e x b/b**2) (x 1.e6 for m/sec)
1873 : ! ui = zonal, vi = meridional, wi = vertical
1874 : !
1875 284544 : do i=lon0,lon1
1876 34429824 : do k=lev0,lev1
1877 34145280 : ui(k,i,j) = -(eey(k,i,j)*zb(i,j)+eez(k,i,j)*xb(i,j))*1.e6_r8/(bmod(i,j)**2)
1878 34145280 : vi(k,i,j) = (eez(k,i,j)*yb(i,j)+eex(k,i,j)*zb(i,j))*1.e6_r8/(bmod(i,j)**2)
1879 34407936 : wi(k,i,j) = (eex(k,i,j)*xb(i,j)-eey(k,i,j)*yb(i,j))*1.e6_r8/(bmod(i,j)**2)
1880 : enddo
1881 : enddo
1882 : !
1883 : ! Output ion drifts in cm/s for oplus_xport call from dpie_coupling:
1884 291840 : do i=lon0,lon1
1885 34407936 : ui(:,i,j) = ui(:,i,j)*100._r8
1886 34407936 : vi(:,i,j) = vi(:,i,j)*100._r8
1887 34429824 : wi(:,i,j) = wi(:,i,j)*100._r8
1888 : enddo
1889 : enddo ! j=lat0,lat1
1890 :
1891 34437120 : call savefld_waccm(eex*100._r8,'EX',nlev,lon0,lon1,lat0,lat1) ! V/m
1892 34437120 : call savefld_waccm(eey*100._r8,'EY',nlev,lon0,lon1,lat0,lat1)
1893 34437120 : call savefld_waccm(eez*100._r8,'EZ',nlev,lon0,lon1,lat0,lat1)
1894 :
1895 7296 : if (debug.and.masterproc) then
1896 : write(iulog,"('ionvel: ion drifts on geo grid: ui=',2e12.4,' vi=',2e12.4,' wi=',2e12.4)") &
1897 0 : minval(ui),maxval(ui), minval(vi),maxval(vi), minval(wi),maxval(wi)
1898 : endif
1899 :
1900 7296 : if (debug_hist) then
1901 0 : if (hist_fld_active('RJAC11')) then
1902 0 : do i=1,nlev
1903 0 : rjac_out(i,lon0:lon1,lat0:lat1) = rjac(lon0:lon1,lat0:lat1,1,1)
1904 : end do
1905 0 : call savefld_waccm(rjac_out,'RJAC11',nlev,lon0,lon1,lat0,lat1)
1906 : endif
1907 :
1908 0 : if (hist_fld_active('RJAC12')) then
1909 0 : do i=1,nlev
1910 0 : rjac_out(i,lon0:lon1,lat0:lat1) = rjac(lon0:lon1,lat0:lat1,1,2)
1911 : end do
1912 0 : call savefld_waccm(rjac_out,'RJAC12',nlev,lon0,lon1,lat0,lat1)
1913 : endif
1914 :
1915 0 : if (hist_fld_active('RJAC21')) then
1916 0 : do i=1,nlev
1917 0 : rjac_out(i,lon0:lon1,lat0:lat1) = rjac(lon0:lon1,lat0:lat1,2,1)
1918 : end do
1919 0 : call savefld_waccm(rjac_out,'RJAC21',nlev,lon0,lon1,lat0,lat1)
1920 : endif
1921 :
1922 0 : if (hist_fld_active('RJAC22')) then
1923 0 : do i=1,nlev
1924 0 : rjac_out(i,lon0:lon1,lat0:lat1) = rjac(lon0:lon1,lat0:lat1,2,2)
1925 : end do
1926 0 : call savefld_waccm(rjac_out,'RJAC22',nlev,lon0,lon1,lat0,lat1)
1927 : endif
1928 : endif
1929 :
1930 7296 : end subroutine ionvel
1931 :
1932 : !-----------------------------------------------------------------------
1933 : end module edynamo
|
