Line data Source code
1 : module edyn_solve
2 : !
3 : ! Prepare stencils and call mudpack PDE solver. This is executed
4 : ! by the root task only, following the gather_edyn call in edynamo.F90.
5 : !
6 : use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals
7 : use cam_logfile ,only: iulog
8 : use edyn_params ,only: finit
9 : use edyn_maggrid ,only: nmlon,nmlonp1,nmlat,nmlath
10 : use edyn_maggrid ,only: res_nlev, res_ngrid
11 : use spmd_utils, only: masterproc
12 : use edyn_solver_coefs, only: nc, cee, cofum
13 :
14 : implicit none
15 : private
16 :
17 : public :: edyn_solve_init
18 : public :: solve_edyn
19 :
20 : !
21 : ! Global 2d fields for root task to complete serial part of dynamo.
22 : ! The zigmxxx, rhs and rims are gathered from subdomains by in sub
23 : ! gather_edyn (edynamo.F90).
24 : !
25 : real(r8),allocatable, dimension(:,:), public :: &
26 : zigm11_glb ,&
27 : zigm22_glb ,&
28 : zigmc_glb ,&
29 : zigm2_glb ,&
30 : rhs_glb
31 : real(r8),allocatable, dimension(:,:,:), public :: &
32 : rim_glb ! pde solver output
33 : real(r8),allocatable, dimension(:,:) :: &
34 : phisolv
35 : !
36 : ! Dimensions of the grid resolutions for the multi-grid PDE:
37 : integer, public, protected :: &
38 : nmlon0, &
39 : nmlat0, &
40 : nmlon1, &
41 : nmlat1, &
42 : nmlon2, &
43 : nmlat2, &
44 : nmlon3, &
45 : nmlat3, &
46 : nmlon4, &
47 : nmlat4, &
48 : nmlon5, &
49 : nmlat5, &
50 : nmlon6, &
51 : nmlat6, &
52 : nmlon7, &
53 : nmlat7
54 : !
55 : ! Space needed for descretized coefficients of of dynamo pde at all levels:
56 : !
57 : integer :: ncee
58 : !
59 : ! The following parameters nc0,nc1,... are pointers to the beginning of
60 : ! the coefficients for each level of resolution.
61 : !
62 : integer :: &
63 : nc0, &
64 : nc1, &
65 : nc2, &
66 : nc3, &
67 : nc4, &
68 : nc5, &
69 : nc6, &
70 : nc7
71 :
72 : real(r8), private, pointer :: &
73 : c0(:), &
74 : c1(:), &
75 : c2(:), &
76 : c3(:), &
77 : c4(:), &
78 : c5(:), &
79 : c6(:), &
80 : c7(:)
81 :
82 : ! phihm is high-latitude potential, set by the high-latitude potential model (e.g. Heelis)
83 : ! or is prescribed (e.g. AMIE)
84 : !
85 : real(r8), allocatable, public :: phihm(:,:) ! high-latitude potential
86 : real(r8), allocatable, public :: pfrac(:,:) ! NH fraction of potential
87 :
88 : contains
89 :
90 : !-----------------------------------------------------------------------
91 768 : subroutine edyn_solve_init
92 : use infnan, only: nan, assignment(=)
93 :
94 3072 : allocate(zigm11_glb(nmlonp1,nmlat))
95 2304 : allocate(zigm22_glb(nmlonp1,nmlat))
96 2304 : allocate(zigmc_glb(nmlonp1,nmlat))
97 2304 : allocate(zigm2_glb(nmlonp1,nmlat))
98 2304 : allocate(rhs_glb(nmlonp1,nmlat))
99 3840 : allocate(rim_glb(nmlonp1,nmlat,2))
100 3072 : allocate(phisolv(0:nmlonp1,0:nmlat+1))
101 :
102 6311424 : phisolv(:,:) = 0._r8
103 :
104 768 : nmlon0=nmlon+1
105 768 : nmlat0=(nmlat +1)/2
106 768 : nmlon1=(nmlon0+1)/2
107 768 : nmlat1=(nmlat0+1)/2
108 768 : nmlon2=(nmlon1+1)/2
109 768 : nmlat2=(nmlat1+1)/2
110 768 : nmlon3=(nmlon2+1)/2
111 768 : nmlat3=(nmlat2+1)/2
112 768 : nmlon4=(nmlon3+1)/2
113 768 : nmlat4=(nmlat3+1)/2
114 768 : nmlon5=(nmlon4+1)/2
115 768 : nmlat5=(nmlat4+1)/2
116 768 : nmlon6=(nmlon5+1)/2
117 768 : nmlat6=(nmlat5+1)/2
118 768 : nmlon7=(nmlon6+1)/2
119 768 : nmlat7=(nmlat6+1)/2
120 :
121 3840 : allocate(cofum(nmlon0,nmlat0,9))
122 :
123 : ncee=10*nmlon0*nmlat0+9*(nmlon1*nmlat1+nmlon2*nmlat2+nmlon3* &
124 768 : nmlat3+nmlon4*nmlat4+nmlon5*nmlat5+nmlon6*nmlat6+nmlon7*nmlat7)
125 :
126 2304 : allocate(cee(ncee))
127 :
128 768 : nc0=1
129 768 : nc1=nc0+10*nmlon0*nmlat0
130 768 : nc2=nc1+9 *nmlon1*nmlat1
131 768 : nc3=nc2+9 *nmlon2*nmlat2
132 768 : nc4=nc3+9 *nmlon3*nmlat3
133 768 : nc5=nc4+9 *nmlon4*nmlat4
134 768 : nc6=nc5+9 *nmlon5*nmlat5
135 768 : nc7=nc6+9 *nmlon6*nmlat6
136 :
137 768 : c0 => cee
138 768 : c1 => cee(nc1:)
139 768 : c2 => cee(nc2:)
140 768 : c3 => cee(nc3:)
141 768 : c4 => cee(nc4:)
142 768 : c5 => cee(nc5:)
143 768 : c6 => cee(nc6:)
144 768 : c7 => cee(nc7:)
145 :
146 2304 : allocate(phihm(nmlonp1,nmlat))
147 3072 : allocate(pfrac(nmlonp1,nmlat0))
148 :
149 768 : phihm = nan
150 768 : pfrac = nan
151 :
152 768 : end subroutine edyn_solve_init
153 :
154 : !-----------------------------------------------------------------------
155 19 : subroutine solve_edyn
156 : !
157 : ! Set up stencils for solver:
158 : !
159 19 : call stencils
160 : !
161 : ! Call mudpack PDE solver:
162 : !
163 19 : call solver(cofum,c0)
164 :
165 19 : end subroutine solve_edyn
166 : !-----------------------------------------------------------------------
167 19 : subroutine stencils
168 : use edyn_params ,only: pi_dyn
169 : use edyn_maggrid,only: dlatm,dlonm
170 : !
171 : ! Locals:
172 : integer :: i,j,jj,jjj,j0,n,ncc,nmaglon,nmaglat, ndx1,ndx2
173 : real(r8) :: sym
174 38 : real(r8) :: cs(nmlat0)
175 :
176 : !
177 : ! Set index array nc and magnetic latitude cosine array:
178 : ! nc pointes to the start of the coefficient array for each level
179 19 : nc(1) = nc0
180 19 : nc(2) = nc1
181 19 : nc(3) = nc2
182 19 : nc(4) = nc3
183 19 : nc(5) = nc4
184 19 : nc(6) = nc5
185 19 : nc(7) = nc6
186 19 : nc(8) = nc7
187 19 : nc(9) = ncee
188 :
189 950 : do j=1,nmlat0
190 950 : cs(j) = cos(pi_dyn/2._r8-(nmlat0-j)*dlatm)
191 : enddo ! j=1,nmlat0
192 : !
193 : ! Set up difference coefficients. Replace zigm11 by A, zigm22 by B,
194 : ! zigmc by C, and zigm2 by D.
195 : !
196 19 : j0 = nmlat0-nmlath
197 950 : do j=1,nmlath ! 1,49 (assuming nmlat=97)
198 931 : jj = nmlath+j-1 ! 49,97
199 931 : jjj = nmlath-j+1 ! 49,1
200 : !
201 : ! factor 4 from 5-point diff. stencil
202 : ! Sigma_(phi lam)/( 4*Delta lam* Delta lon )
203 : ! Sigma_(phi lam)/( 4*Delta lam* Delta lon )
204 : ! Sigma_(lam lam)*cos(lam_m)*DT0DTS/(Delta lam)^2
205 : ! -zigmc_north = southern hemis. 49,1 equator-pole
206 : ! -zigm2_north = southern hemis. 49,1 equator-pole
207 : ! zigm22 = southern hemis. 49,1 equator-pole
208 : !
209 76342 : do i=1,nmlonp1
210 0 : zigmc_glb(i,jj) = (zigmc_glb(i,jj)+zigm2_glb(i,jj))/ &
211 75411 : (4._r8*dlatm*dlonm)
212 0 : zigm2_glb(i,jj) = zigmc_glb(i,jj)-2._r8*zigm2_glb(i,jj)/ &
213 75411 : (4._r8*dlatm*dlonm)
214 75411 : zigm22_glb(i,jj) = zigm22_glb(i,jj)*cs(j0+j)/dlatm**2
215 75411 : zigmc_glb(i,jjj) = -zigmc_glb(i,jj)
216 75411 : zigm2_glb(i,jjj) = -zigm2_glb(i,jj)
217 76342 : zigm22_glb(i,jjj) = zigm22_glb(i,jj)
218 : enddo ! i=1,nmlonp1
219 950 : if (j /= nmlath) then
220 : !
221 : ! Sigma_(phi phi)/( cos(lam_m)*DT0DTS*(Delta lon)^2 )
222 : ! zigm11 = southern hemis. 49,1 equator-pole
223 : !
224 74784 : do i = 1,nmlonp1
225 73872 : zigm11_glb(i,jj) = zigm11_glb(i,jj)/(cs(j0+j)*dlonm**2)
226 74784 : zigm11_glb(i,jjj) = zigm11_glb(i,jj)
227 : enddo
228 : endif
229 : enddo ! j=1,nmlath
230 : !
231 : ! Set zigm11 to zero at megnetic poles to avoid floating exception
232 : ! (values at poles are not used):
233 : !
234 1558 : do i = 1,nmlonp1
235 1539 : zigm11_glb(i,1) = 0.0_r8
236 1558 : zigm11_glb(i,nmlat) = 0.0_r8
237 : enddo
238 :
239 : ! Clear array for difference stencils at all levels:
240 19 : call clearcee(cee,nmlon0,nmlat0)
241 : !
242 : ! Init cofum coefficients:
243 687268 : cofum(:,:,:) = finit
244 : !
245 : ! Calculate contribution to stencils from each PDE coefficient
246 : !
247 : ! Sigma_(phi phi)/( cos(lam_m)*dt0dts*(Delta lon)^2 )
248 19 : sym = 1._r8
249 19 : call stencmd(zigm11_glb,nmlon0,nmlat0,sym,cee,1)
250 : !
251 : ! Sigma_(lam lam)*cos(lam_m)*dt0dts/(Delta lam)^2
252 19 : sym = 1._r8
253 19 : call stencmd(zigm22_glb,nmlon0,nmlat0,sym,cee,4)
254 : !
255 : ! Sigma_(phi lam)/( 4*Delta lam* Delta lon )
256 19 : sym = -1._r8
257 19 : call stencmd(zigmc_glb,nmlon0,nmlat0,sym,cee,2)
258 : !
259 : ! Sigma_(lam phi)/( 4*Delta lam* Delta lon )
260 19 : sym = -1._r8
261 19 : call stencmd(zigm2_glb,nmlon0,nmlat0,sym,cee,3)
262 : !
263 : ! Insert RHS in finest stencil:
264 950 : do j = 1,nmlat0
265 931 : jj = nmlath-nmlat0+j
266 76361 : do i = 1,nmlon0
267 75411 : ndx1 = 9*nmlat0*nmlon0 + (j-1)*nmlon0 + i
268 76342 : c0(ndx1) = rhs_glb(i,jj)
269 : enddo ! i = 1,nmlon0
270 : enddo ! j = 1,nmlat0
271 19 : ndx1 = 9*nmlat0*nmlon0 + nmlonp1
272 19 : ndx2 = 9*nmlat0*nmlon0 + 1
273 19 : c0(ndx1) = c0(ndx2)
274 : !
275 : ! Set boundary condition at the pole:
276 19 : call edges(c0,nmlon0,nmlat0)
277 19 : call edges(c1,nmlon1,nmlat1)
278 19 : call edges(c2,nmlon2,nmlat2)
279 19 : call edges(c3,nmlon3,nmlat3)
280 19 : call edges(c4,nmlon4,nmlat4)
281 19 : if ( res_nlev > 5 ) then
282 0 : call edges(c5,nmlon5,nmlat5)
283 : endif
284 19 : if ( res_nlev > 6 ) then
285 0 : call edges(c6,nmlon6,nmlat6)
286 : endif
287 19 : if ( res_nlev > 7 ) then
288 0 : call edges(c7,nmlon7,nmlat7)
289 : endif
290 19 : call edges(cofum,nmlon0,nmlat0)
291 : !
292 : ! Divide stencils by cos(lam_0) (not rhs):
293 19 : call divide(c0,nmlon0,nmlat0,nmlon0,cs,1)
294 19 : call divide(c1,nmlon1,nmlat1,nmlon0,cs,1)
295 19 : call divide(c2,nmlon2,nmlat2,nmlon0,cs,1)
296 19 : call divide(c3,nmlon3,nmlat3,nmlon0,cs,1)
297 19 : call divide(c4,nmlon4,nmlat4,nmlon0,cs,1)
298 19 : if ( res_nlev > 5 ) then
299 0 : call divide(c5,nmlon5,nmlat5,nmlon0,cs,1)
300 : endif
301 19 : if ( res_nlev > 6 ) then
302 0 : call divide(c6,nmlon6,nmlat6,nmlon0,cs,1)
303 : endif
304 19 : if ( res_nlev > 7 ) then
305 0 : call divide(c7,nmlon7,nmlat7,nmlon0,cs,1)
306 : endif
307 19 : call divide(cofum,nmlon0,nmlat0,nmlon0,cs,0)
308 : !
309 : ! Set value of solution to 1. at pole:
310 1558 : do i=1,nmlon0
311 1539 : ndx1 = 9*nmlat0*nmlon0 + (nmlat0-1)*nmlon0 + i
312 1558 : c0(ndx1) = 1._r8
313 : enddo
314 : !
315 : ! Modify stencils and RHS so that the NH high lat potential is inserted at
316 : ! high latitude. The SH high lat potential will be added back later.
317 : ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat
318 : ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg,
319 : ! dynamo < 60 deg, and combination between 60-75 mag lat.
320 : ! The dynamo is symmetric about the magnetic equator, but the high latitude
321 : ! is anti-symmetric in both hemispheres. However, since Mudpack uses the
322 : ! NH potential pattern, then the SH potential pattern must be added
323 : ! back into the 2-D phim before the call threed, and before it is
324 : ! transformed to geographic coordinates.
325 : !
326 19 : ncc = 1
327 19 : nmaglon = nmlon0
328 19 : nmaglat = nmlat0
329 114 : do n=1,res_nlev ! resolution levels
330 95 : call stenmd(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac)
331 95 : ncc = ncc+9*nmaglon*nmaglat
332 95 : if (n==1) ncc = ncc+nmaglon*nmaglat ! rhs is in 10th slot
333 95 : nmaglon = (nmaglon+1)/2
334 209 : nmaglat = (nmaglat+1)/2
335 : enddo
336 :
337 19 : end subroutine stencils
338 : !-----------------------------------------------------------------------
339 19 : subroutine clearcee(cee,nlon0,nlat0)
340 : !
341 : ! Zero C arrays for stencil coefficients.
342 : ! Cee will contain:
343 : ! c0(nmlon0,nmlat0,10), c1(nmlon1,nmlat1,9), c2(nmlon2,nmlat2,9),
344 : ! c3(nmlon3,nmlat3,9), c4(nmlon4,nmlat4,9)
345 : !
346 : ! Args:
347 : integer,intent(in) :: nlon0,nlat0
348 : real(r8),intent(out) :: cee(*)
349 : !
350 : ! Local:
351 : integer :: nlon,nlat,n,m,i
352 : !
353 : ! Compute total size of cee
354 19 : nlon = nlon0
355 19 : nlat = nlat0
356 19 : n = 0
357 114 : do m=1,res_nlev ! resolution levels
358 95 : n = n+nlon*nlat
359 95 : nlon = (nlon+1)/2
360 114 : nlat = (nlat+1)/2
361 : enddo
362 19 : n = 9*n+nlon0*nlat0
363 : !
364 : ! Clear cee:
365 993358 : do i=1,n
366 993358 : cee(i) = 0._r8
367 : enddo
368 19 : end subroutine clearcee
369 : !-----------------------------------------------------------------------
370 76 : subroutine stencmd(zigm,nlon0,nlat0,sym,cee,ncoef)
371 : !
372 : ! Calculate contribution fo 3 by 3 stencil from coefficient zigm
373 : ! at each grid point and level.
374 : !
375 : ! Args:
376 : integer,intent(in) :: &
377 : nlon0, & ! longitude dimension of finest grid level
378 : nlat0, & ! latitude dimension of finest grid level
379 : ncoef ! integer identifier of coefficient
380 : real(r8),intent(in) :: &
381 : zigm(nlon0,nlat0), & ! coefficients (nlon0+1/2,(nlat0+1)/2)
382 : sym ! 1. if zigm symmetric w.r.t. equator, -1 otherwise
383 : real(r8),intent(inout) :: & ! output stencil array consisting of c0,c1,c2,c3,c4
384 : cee(*)
385 : !
386 : ! Local:
387 : integer :: nc,nlon,nlat,n
388 152 : real(r8) :: wkarray(-res_ngrid+1:nmlon0+res_ngrid,nmlat0)
389 : !
390 : ! Perform half-way interpolation and extend zigm in wkarray:
391 : !
392 76 : call htrpex(zigm,nlon0,nlat0,sym,wkarray)
393 : !
394 : ! Calculate contribution to stencil for each grid point and level:
395 : !
396 76 : nc = 1
397 76 : nlon = nlon0
398 76 : nlat = nlat0
399 : !
400 : ! Calculate modified and unmodified stencil on finest grid
401 : !
402 76 : call cnmmod(nlon0,nlon,nlat,cee(nc),ncoef,wkarray,cofum)
403 : !
404 : ! Stencils on other grid levels remain the same.
405 76 : nc = nc+10*nlon*nlat
406 76 : nlon = (nlon+1)/2
407 76 : nlat = (nlat+1)/2
408 : !
409 380 : do n=2,res_nlev
410 304 : call cnm(nlon0,nlon,nlat,cee(nc),ncoef,wkarray)
411 304 : nc = nc+9*nlon*nlat
412 304 : if (n==1) nc = nc+nlon*nlat
413 304 : nlon = (nlon+1)/2
414 684 : nlat = (nlat+1)/2
415 : enddo ! n=1,5
416 76 : end subroutine stencmd
417 : !-----------------------------------------------------------------------
418 76 : subroutine htrpex(coeff,nmlon0,nmlat0,sym,wkarray)
419 : !
420 : ! Perform half-way interpolation on array coeff and extend over 16 grid
421 : ! points. Result returned in wkarray.
422 : !
423 : ! Args:
424 : integer,intent(in) :: nmlon0,nmlat0
425 : real(r8),intent(in) :: coeff(nmlon0,nmlat0),sym
426 : real(r8),intent(out) :: wkarray(-res_ngrid+1:nmlon0+res_ngrid,nmlat0)
427 : !
428 : ! Local:
429 : integer :: i,j,jj
430 : !
431 : ! Copy coeff into positions in wkarray:
432 3800 : do j=1,nmlat0
433 3724 : jj = nmlat0-j+1
434 305444 : do i=1,nmlon0
435 305368 : wkarray(i,j) = sym*coeff(i,jj)
436 : enddo ! i=1,nmlon0
437 : enddo ! j=1,nmlat0
438 : !
439 : ! Extend over 2*res_ngrid grid spaces to allow for a total of res_nlev grid levels:
440 1292 : do i=1,res_ngrid
441 60876 : do j=1,nmlat0
442 59584 : wkarray(1-i,j) = wkarray(nmlon0-i,j)
443 60800 : wkarray(nmlon0+i,j) = wkarray(1+i,j)
444 : enddo ! j=1,nmlat0
445 : enddo ! i=1,res_ngrid
446 76 : end subroutine htrpex
447 : !-----------------------------------------------------------------------
448 304 : subroutine cnm(nlon0,nlon,nlat,c,ncoef,wkarray)
449 : !
450 : ! Compute contribution to stencil from zigm(ncoef) on grid nlon by nlat,
451 : ! Finest grid is nlon0.
452 : !
453 : ! Args:
454 : integer,intent(in) :: &
455 : nlon0, & ! finest grid dimensions
456 : nlon,nlat ! output grid dimensions
457 : real(r8),intent(in) :: wkarray(-res_ngrid+1:nmlon0+res_ngrid,nmlat0)
458 : !
459 : ! ncoef: integer id of coefficient:
460 : ! ncoef = 1 for zigm11
461 : ! ncoef = 2 for zigm12 (=zigmc+zigm2)
462 : ! ncoef = 3 for zigm21 (=zigmc-zigm2)
463 : ! ncoef = 4 for zigm22
464 : !
465 : integer,intent(in) :: ncoef
466 : real(r8),intent(inout) :: &
467 : c(nlon,nlat,*) ! output array for grid point stencils at resolution nlon x nlat
468 : !
469 : ! Local:
470 : integer :: i,j,nint,i0,j0
471 : ! For now, retain this pi to insure bit compatability w/ old code
472 : real(r8),parameter :: pi=3.141592654_r8
473 608 : real(r8) :: wk(nlon0,3)
474 : !
475 : ! Compute separation of grid points of resolution nlon x nlat within
476 : ! grid of resolution nlon0. Evaluate dlon and dlat, grid spacing
477 : ! of nlon x nlat.
478 : !
479 304 : nint = (nlon0-1)/(nlon-1)
480 : !
481 : ! Scan wkarray nlon x nlat calculating and adding contributions to stencil
482 : ! from zigm(ncoef)
483 304 : i0 = 1-nint
484 304 : j0 = 1-nint
485 : !
486 : ! zigm11:
487 : ! am 2001-6-27 include boundary condition at equator
488 304 : if (ncoef==1) then
489 931 : do j = 1,nlat-1
490 26011 : do i = 1,nlon
491 100320 : c(i,j,1) = c(i,j,1)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
492 125400 : wkarray(i0+(i+1)*nint,j0+j*nint))
493 : c(i,j,5) = c(i,j,5)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
494 25080 : wkarray(i0+(i-1)*nint,j0+j*nint))
495 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+(i+1)*nint,j0+j*nint)+ &
496 25935 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+wkarray(i0+(i-1)*nint,j0+j*nint))
497 : enddo ! i = 1,nlon
498 : enddo ! j = 2,nlat-1
499 : !
500 : ! zigm12 (=zigmc+zigm2)
501 228 : elseif (ncoef==2) then
502 855 : do j = 2,nlat-1
503 24434 : do i = 1,nlon
504 94316 : c(i,j,2) = c(i,j,2)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
505 117895 : wkarray(i0+(i+1)*nint,j0+j*nint))
506 : c(i,j,4) = c(i,j,4)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
507 23579 : wkarray(i0+(i-1)*nint,j0+j*nint))
508 : c(i,j,6) = c(i,j,6)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
509 23579 : wkarray(i0+(i-1)*nint,j0+j*nint))
510 : c(i,j,8) = c(i,j,8)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
511 23579 : wkarray(i0+(i+1)*nint,j0+j*nint))
512 23579 : wk(i,1) = 0.5_r8*(wkarray(i0+(i+1)*nint,j0+j*nint)- &
513 23579 : wkarray(i0+(i-1)*nint,j0+j*nint))
514 23579 : wk(i,2) = (c(i,j,3)+wk(i,1))*(c(i,j,7)-wk(i,1))
515 23579 : wk(i,3) = sign(wk(i,1),c(i,j,3)+c(i,j,7))
516 23579 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
517 23579 : c(i,j,3) = c(i,j,3)+wk(i,1)+wk(i,3)
518 23579 : c(i,j,7) = c(i,j,7)-wk(i,1)+wk(i,3)
519 24358 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
520 : enddo ! i = 1,nlon
521 : enddo ! j = 2,nlat-1
522 : !
523 : ! zigm21 (=zigmc-zigm2)
524 152 : elseif (ncoef==3) then
525 855 : do j = 2,nlat-1
526 24434 : do i = 1,nlon
527 94316 : c(i,j,2) = c(i,j,2)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
528 117895 : wkarray(i0+i*nint,j0+(j+1)*nint))
529 : c(i,j,4) = c(i,j,4)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
530 23579 : wkarray(i0+i*nint,j0+(j+1)*nint))
531 : c(i,j,6) = c(i,j,6)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
532 23579 : wkarray(i0+i*nint,j0+(j-1)*nint))
533 : c(i,j,8) = c(i,j,8)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
534 23579 : wkarray(i0+i*nint,j0+(j-1)*nint))
535 23579 : wk(i,1) = 0.5_r8*(wkarray(i0+i*nint,j0+(j+1)*nint)- &
536 23579 : wkarray(i0+i*nint,j0+(j-1)*nint))
537 23579 : wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1))
538 23579 : wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5))
539 23579 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
540 23579 : c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3)
541 23579 : c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3)
542 24358 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
543 : enddo ! i = 1,nlon
544 : enddo ! j = 2,nlat-1
545 : !
546 : ! Low latitude boundary condition:
547 1577 : j = 1
548 1577 : do i=1,nlon
549 4503 : c(i,j,2) = c(i,j,2)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
550 4503 : wkarray(i0+i*nint,j0+(j+1)*nint))
551 : c(i,j,4) = c(i,j,4)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
552 1501 : wkarray(i0+i*nint,j0+(j+1)*nint))
553 1501 : wk(i,1) = 0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
554 1501 : wkarray(i0+i*nint,j0+(j+1)*nint))
555 1501 : wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1))
556 1501 : wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5))
557 1501 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
558 1501 : c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3)
559 1501 : c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3)
560 1577 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
561 : enddo ! i=1,nlon
562 : !
563 : ! zigm22:
564 76 : elseif (ncoef==4) then
565 855 : do j = 2,nlat-1
566 24434 : do i = 1,nlon
567 94316 : c(i,j,3) = c(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
568 117895 : wkarray(i0+i*nint,j0+(j+1)*nint))
569 : c(i,j,7) = c(i,j,7)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
570 23579 : wkarray(i0+i*nint,j0+(j-1)*nint))
571 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+(j-1)*nint)+ &
572 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+ &
573 24358 : wkarray(i0+i*nint,j0+(j+1)*nint))
574 : enddo ! i = 1,nlon
575 : enddo ! j = 2,nlat-1
576 : !
577 : ! Low latitude boundary condition:
578 1577 : j = 1
579 1577 : do i=1,nlon
580 4503 : c(i,j,3) = c(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
581 4503 : wkarray(i0+i*nint,j0+(j+1)*nint))
582 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
583 1577 : wkarray(i0+i*nint,j0+(j+1)*nint))
584 : enddo ! i=1,nlon
585 : endif ! ncoef
586 304 : end subroutine cnm
587 : !-----------------------------------------------------------------------
588 76 : subroutine cnmmod(nlon0,nlon,nlat,c,ncoef,wkarray,cofum)
589 : !
590 : ! Compute contribution to stencil from zigm(ncoef) on grid nlon by nlat,
591 : ! Finest grid is nlon0.
592 : !
593 : ! Args:
594 : integer,intent(in) :: &
595 : nlon0, & ! finest grid dimensions
596 : nlon,nlat ! output grid dimensions
597 : real(r8),intent(in) :: wkarray(-res_ngrid+1:nmlon0+res_ngrid,nmlat0)
598 : real(r8),dimension(nmlon0,nmlat0,9),intent(inout) :: cofum
599 : !
600 : ! ncoef: integer id of coefficient:
601 : ! ncoef = 1 for zigm11
602 : ! ncoef = 2 for zigm12 (=zigmc+zigm2)
603 : ! ncoef = 3 for zigm21 (=zigmc-zigm2)
604 : ! ncoef = 4 for zigm22
605 : !
606 : integer,intent(in) :: ncoef
607 : real(r8),intent(inout) :: &
608 : c(nlon,nlat,*) ! output array for grid point stencils at resolution nlon x nlat
609 : !
610 : ! Local:
611 : integer :: i,j,nint,i0,j0
612 : ! For now, retain this pi to insure bit compatability w/ old code
613 : real(r8),parameter :: pi=3.141592654_r8
614 152 : real(r8) :: wk(nlon0,3)
615 : !
616 : ! Compute separation of grid points of resolution nlon x nlat within
617 : ! grid of resolution nlon0. Evaluate dlon and dlat, grid spacing
618 : ! of nlon x nlat.
619 : !
620 76 : nint = (nlon0-1)/(nlon-1)
621 : !
622 : ! Scan wkarray nlon x nlat calculating and adding contributions to stencil
623 : ! from zigm(ncoef)
624 76 : i0 = 1-nint
625 76 : j0 = 1-nint
626 : !
627 : ! zigm11:
628 : ! am 2001-6-27 include boundary condition at equator
629 76 : if (ncoef==1) then
630 931 : do j = 1,nlat-1
631 74803 : do i = 1,nlon
632 295488 : c(i,j,1) = c(i,j,1)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
633 369360 : wkarray(i0+(i+1)*nint,j0+j*nint))
634 : c(i,j,5) = c(i,j,5)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
635 73872 : wkarray(i0+(i-1)*nint,j0+j*nint))
636 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+(i+1)*nint,j0+j*nint)+ &
637 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+ &
638 73872 : wkarray(i0+(i-1)*nint,j0+j*nint))
639 : !
640 : ! Unmodified:
641 73872 : cofum(i,j,1) = cofum(i,j,1)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
642 73872 : wkarray(i0+(i+1)*nint,j0+j*nint))
643 : cofum(i,j,5) = cofum(i,j,5)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
644 73872 : wkarray(i0+(i-1)*nint,j0+j*nint))
645 : cofum(i,j,9) = cofum(i,j,9)-0.5_r8*(wkarray(i0+(i+1)*nint,j0+j*nint)+ &
646 74784 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+wkarray(i0+(i-1)*nint,j0+j*nint))
647 : enddo ! i = 1,nlon
648 : enddo ! j = 2,nlat-1
649 : !
650 : ! zigm12 (=zigmc+zigm2)
651 57 : elseif (ncoef==2) then
652 912 : do j = 2,nlat-1
653 73245 : do i = 1,nlon
654 289332 : c(i,j,2) = c(i,j,2)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
655 361665 : wkarray(i0+(i+1)*nint,j0+j*nint))
656 : c(i,j,4) = c(i,j,4)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
657 72333 : wkarray(i0+(i-1)*nint,j0+j*nint))
658 : c(i,j,6) = c(i,j,6)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
659 72333 : wkarray(i0+(i-1)*nint,j0+j*nint))
660 : c(i,j,8) = c(i,j,8)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
661 72333 : wkarray(i0+(i+1)*nint,j0+j*nint))
662 72333 : wk(i,1) = 0.5_r8*(wkarray(i0+(i+1)*nint,j0+j*nint)- &
663 72333 : wkarray(i0+(i-1)*nint,j0+j*nint))
664 : !
665 : ! Unmodified:
666 72333 : cofum(i,j,2) = c(i,j,2)
667 72333 : cofum(i,j,4) = c(i,j,4)
668 72333 : cofum(i,j,6) = c(i,j,6)
669 72333 : cofum(i,j,8) = c(i,j,8)
670 72333 : cofum(i,j,3) = cofum(i,j,3)+wk(i,1)
671 72333 : cofum(i,j,7) = cofum(i,j,7)-wk(i,1)
672 : !
673 72333 : wk(i,2) = (c(i,j,3)+wk(i,1))*(c(i,j,7)-wk(i,1))
674 72333 : wk(i,3) = sign(wk(i,1),c(i,j,3)+c(i,j,7))
675 72333 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
676 72333 : c(i,j,3) = c(i,j,3)+wk(i,1)+wk(i,3)
677 72333 : c(i,j,7) = c(i,j,7)-wk(i,1)+wk(i,3)
678 73226 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
679 : enddo ! i = 1,nlon
680 : enddo ! j = 2,nlat-1
681 : !
682 : ! zigm21 (=zigmc-zigm2)
683 38 : elseif (ncoef==3) then
684 912 : do j = 2,nlat-1
685 73245 : do i = 1,nlon
686 289332 : c(i,j,2) = c(i,j,2)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
687 361665 : wkarray(i0+i*nint,j0+(j+1)*nint))
688 : c(i,j,4) = c(i,j,4)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
689 72333 : wkarray(i0+i*nint,j0+(j+1)*nint))
690 : c(i,j,6) = c(i,j,6)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
691 72333 : wkarray(i0+i*nint,j0+(j-1)*nint))
692 : c(i,j,8) = c(i,j,8)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
693 72333 : wkarray(i0+i*nint,j0+(j-1)*nint))
694 72333 : wk(i,1) = 0.5_r8*(wkarray(i0+i*nint,j0+(j+1)*nint)- &
695 72333 : wkarray(i0+i*nint,j0+(j-1)*nint))
696 : !
697 : ! Unmodified:
698 72333 : cofum(i,j,2) = c(i,j,2)
699 72333 : cofum(i,j,4) = c(i,j,4)
700 72333 : cofum(i,j,6) = c(i,j,6)
701 72333 : cofum(i,j,8) = c(i,j,8)
702 72333 : cofum(i,j,1) = cofum(i,j,1)+wk(i,1)
703 72333 : cofum(i,j,5) = cofum(i,j,5)-wk(i,1)
704 : !
705 72333 : wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1))
706 72333 : wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5))
707 72333 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
708 72333 : c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3)
709 72333 : c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3)
710 73226 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
711 : enddo ! i = 1,nlon
712 : enddo ! j = 2,nlat-1
713 : !
714 : ! Low latitude boundary condition:
715 1558 : j = 1
716 1558 : do i=1,nlon
717 4617 : c(i,j,2) = c(i,j,2)+.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
718 4617 : wkarray(i0+i*nint,j0+(j+1)*nint))
719 : c(i,j,4) = c(i,j,4)-.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
720 1539 : wkarray(i0+i*nint,j0+(j+1)*nint))
721 1539 : wk(i,1) = .5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
722 1539 : wkarray(i0+i*nint,j0+(j+1)*nint))
723 :
724 1539 : cofum(i,j,2) = c(i,j,2)
725 1539 : cofum(i,j,4) = c(i,j,4)
726 1539 : cofum(i,j,1) = cofum(i,j,1)+wk(i,1)
727 1539 : cofum(i,j,5) = cofum(i,j,5)-wk(i,1)
728 :
729 1539 : wk(i,2) = (c(i,j,1)+wk(i,1))*(c(i,j,5)-wk(i,1))
730 1539 : wk(i,3) = sign(wk(i,1),c(i,j,1)+c(i,j,5))
731 1539 : if (wk(i,2) >= 0._r8) wk(i,3) = 0._r8
732 1539 : c(i,j,1) = c(i,j,1)+wk(i,1)+wk(i,3)
733 1539 : c(i,j,5) = c(i,j,5)-wk(i,1)+wk(i,3)
734 1558 : c(i,j,9) = c(i,j,9)-2._r8*wk(i,3)
735 : enddo ! i=1,nlon
736 : !
737 : ! zigm22:
738 19 : elseif (ncoef==4) then
739 912 : do j = 2,nlat-1
740 73245 : do i = 1,nlon
741 289332 : c(i,j,3) = c(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
742 361665 : wkarray(i0+i*nint,j0+(j+1)*nint))
743 : c(i,j,7) = c(i,j,7)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
744 72333 : wkarray(i0+i*nint,j0+(j-1)*nint))
745 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+(j-1)*nint)+ &
746 72333 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+wkarray(i0+i*nint,j0+(j+1)*nint))
747 : !
748 : ! Unmodified:
749 72333 : cofum(i,j,3) = cofum(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
750 72333 : wkarray(i0+i*nint,j0+(j+1)*nint))
751 : cofum(i,j,7) = cofum(i,j,7)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
752 72333 : wkarray(i0+i*nint,j0+(j-1)*nint))
753 : cofum(i,j,9) = cofum(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+(j-1)*nint)+ &
754 73226 : 2._r8*wkarray(i0+i*nint,j0+j*nint)+wkarray(i0+i*nint,j0+(j+1)*nint))
755 : enddo ! i = 1,nlon
756 : enddo ! j = 2,nlat-1
757 : !
758 : ! Low latitude boundary condition:
759 1558 : j = 1
760 1558 : do i=1,nlon
761 4617 : c(i,j,3) = c(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
762 4617 : wkarray(i0+i*nint,j0+(j+1)*nint))
763 : c(i,j,9) = c(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
764 1539 : wkarray(i0+i*nint,j0+(j+1)*nint))
765 1539 : cofum(i,j,3) = cofum(i,j,3)+0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
766 1539 : wkarray(i0+i*nint,j0+(j+1)*nint))
767 : cofum(i,j,9) = cofum(i,j,9)-0.5_r8*(wkarray(i0+i*nint,j0+j*nint)+ &
768 1558 : wkarray(i0+i*nint,j0+(j+1)*nint))
769 : enddo ! i=1,nlon
770 : endif ! ncoef
771 76 : end subroutine cnmmod
772 : !--------------------------------------------------------------------
773 114 : subroutine edges(c,nlon,nlat)
774 : !
775 : ! Insert equatorial and polar boundary conditions in stencil c(nlon,nlat,9)
776 : !
777 : ! Args:
778 : integer,intent(in) :: nlon,nlat
779 : real(r8),intent(out) :: c(*)
780 : !
781 : ! Local:
782 : integer :: n,i, ndx
783 :
784 1026 : do n=1,8
785 37658 : do i=1,nlon
786 36632 : ndx = (n-1)*nlat*nlon + (nlat-1)*nlon + i
787 37544 : c(ndx) = 0._r8
788 : enddo
789 : enddo
790 4693 : do i=1,nlon
791 4579 : ndx = 8*nlat*nlon + (nlat-1)*nlon + i
792 4693 : c(ndx) = 1._r8
793 : enddo
794 114 : end subroutine edges
795 : !--------------------------------------------------------------------
796 114 : subroutine divide(c,nlon,nlat,nlon0,cs,igrid)
797 : !
798 : ! Divide stencil C by cos(theta(i,j))
799 : !
800 : ! Args:
801 : integer,intent(in) :: nlon,nlat,nlon0,igrid
802 : real(r8),intent(in) :: cs(:)
803 : real(r8),intent(out) :: c(*)
804 : !
805 : ! Local:
806 : integer :: nint,j0,n,j,i, ndx
807 : !
808 114 : nint = (nlon0-1)/(nlon-1)
809 114 : j0 = 1-nint
810 1140 : do n = 1,9
811 25251 : do j = 1,nlat-1
812 1580553 : do i = 1,nlon
813 1555416 : ndx = (n-1)*nlat*nlon + (j-1)*nlon + i
814 1579527 : c(ndx) = c(ndx)/(cs(j0+j*nint)*nint**2)
815 : enddo ! i = 1,nlon
816 : enddo ! j = 1,nlat-1
817 : enddo ! n = 1,9
818 : !
819 114 : if (nint==1.and.igrid > 0) then
820 1558 : do i = 1,nlon
821 1539 : ndx = 9*nlat*nlon + i
822 1558 : c(ndx) = c(ndx)/cs(1)
823 : enddo ! i = 1,nlon
824 : endif
825 114 : end subroutine divide
826 : !--------------------------------------------------------------------
827 95 : subroutine stenmd(inlon,inlat,c,phihm,pfrac)
828 : use edyn_params ,only: dtr
829 : use edyn_maggrid,only: dlatm
830 : !
831 : ! Modify stencil to set potential to heelis value within auroral circle.
832 : !
833 : ! Args:
834 : integer,intent(in) :: inlon,inlat
835 : real(r8),intent(inout) :: c(inlon,inlat,*)
836 : real(r8),dimension(nmlon0,nmlat0),intent(in) :: &
837 : phihm, & ! heelis potential (from subs potm, flwv32)
838 : pfrac ! fractional presence of dynamo (from sub colath)
839 : !
840 : ! Local:
841 : integer :: nint,i0,j0,i,j,n,jj
842 : real(r8) :: real8
843 : !
844 : ! Compute separation of grid points for this resolution:
845 95 : nint = (nmlon0-1)/(inlon-1)
846 95 : i0 = 1-nint
847 95 : j0 = 1-nint
848 : !
849 : ! If nint==1, then we are at the highest resolution.
850 : ! Correct RHS, which is in c(10)
851 : !
852 95 : if (nint==1) then
853 950 : do j=1,inlat
854 76361 : do i=1,inlon
855 150822 : c(i,j,10) = pfrac(i,j)*c(i,j,10)+(1._r8-pfrac(i,j))*c(i,j,9)* &
856 227164 : (dlatm/(10._r8*dtr))**2*phihm(i,j)
857 : enddo ! i=1,inlon
858 : enddo ! j=1,inlat
859 : endif
860 : !
861 : ! Modify stencil, c(i,j,n),n=1,9:
862 : !
863 95 : real8 = dble(nint)
864 : if (nint==1) then
865 950 : do j=1,inlat
866 931 : jj = j0+j*nint
867 8379 : do n = 1,8
868 611667 : do i = 1,inlon
869 603288 : c(i,j,n) = c(i,j,n)*pfrac(i0+i*nint,jj)
870 610736 : cofum(i,j,n) = cofum(i,j,n)*pfrac(i0+i*nint,jj)
871 : enddo ! i = 1,inlon
872 : enddo ! n = 1,8
873 76361 : do i = 1,inlon
874 301644 : c(i,j,9) = c(i,j,9)*pfrac(i0+i*nint,jj)+ &
875 : (1._r8-pfrac(i0+i*nint,jj))*c(i,j,9)* &
876 301644 : (dlatm*real8/(10._r8*dtr))**2
877 0 : cofum(i,j,9) =cofum(i,j,9)*pfrac(i0+i*nint,jj)+ &
878 : (1._r8-pfrac(i0+i*nint,jj))*cofum(i,j,9)* &
879 76342 : (dlatm*real8/(10._r8*dtr))**2
880 : enddo ! i = 1,inlon
881 : enddo ! j=1,inlat
882 : else ! nint /= 1
883 1007 : do j=1,inlat
884 931 : jj = j0+j*nint
885 8379 : do n = 1,8
886 221027 : do i = 1,inlon
887 220096 : c(i,j,n) = c(i,j,n)*pfrac(i0+i*nint,jj)
888 : enddo ! i = 1,inlon
889 : enddo ! n = 1,8
890 27588 : do i = 1,inlon
891 106324 : c(i,j,9) = c(i,j,9)*pfrac(i0+i*nint,jj)+ &
892 : (1._r8-pfrac(i0+i*nint,jj))*c(i,j,9)* &
893 133836 : (dlatm*real8/(10._r8*dtr))**2
894 : enddo ! i = 1,inlon
895 : enddo ! j=1,inlat
896 : endif ! nint
897 95 : end subroutine stenmd
898 : !--------------------------------------------------------------------
899 19 : subroutine solver(cofum,c0)
900 : use edyn_mudmod, only: mudmod
901 : use edyn_muh2cr, only: muh
902 : !
903 : ! Call mudpack to solve PDE. Solution is returned in rim:
904 : ! real,dimension(nmlonp1,nmlat,2) :: rim
905 : ! Mudpack solvers:
906 : ! isolve = 0 org. mud v5. (modified stencils neq direct solution)
907 : ! isolve = 1 muh hybrid solver (no convergence => only as direct solver)
908 : ! isolve = 2 modified mud (residual calculated with unmodified stencils
909 : ! same solution as with direct solver, if same
910 : ! coefficient matrix is used)
911 : ! Only isolve==2 is supported in edynamo.
912 : !
913 : ! Args:
914 : real(r8),dimension(nmlon0,nmlat0,9),intent(in) :: cofum
915 : real(r8),intent(in) :: c0(nmlon0,nmlat0,10)
916 : !
917 : ! Local:
918 : integer :: i,j,jntl,ier,isolve
919 : real(r8) :: l2norm
920 38 : real(r8),dimension(nmlon0,nmlat0) :: ressolv
921 38 : real(r8),dimension(nmlon0,nmlat0,9) :: cofum_solv
922 :
923 : ! Module data:
924 : ! real,dimension(nmlonp1,nmlat,2) :: rim_glb ! pde solver output
925 :
926 19 : jntl = 0
927 19 : ier = 0
928 19 : isolve = 2
929 19 : call mudmod(rim_glb,phisolv,jntl,isolve,res_nlev,ier)
930 19 : if (ier < 0 ) then ! not converged
931 0 : if (masterproc) write(iulog,*) 'solver: use muh direct solver'
932 0 : call muh(rim_glb,nmlon,nmlat,res_nlev,jntl)
933 : endif
934 :
935 19 : l2norm=0._r8
936 76361 : ressolv = 0.0_r8
937 950 : do j = 1,nmlat0
938 75430 : do i = 1,nmlon0-1
939 744800 : cofum_solv(i,j,:) = cofum(i,j,:)
940 : !
941 : ! fields: phisolv(0:nmlonp1,0:nmlat+1) ! 2d solution/ electric potential
942 : !
943 : ressolv(i,j) = ( &
944 0 : cofum_solv(i,j,1)*phisolv(i+1,j)+ &
945 74480 : cofum_solv(i,j,2)*phisolv(i+1,j+1)+ &
946 74480 : cofum_solv(i,j,3)*phisolv(i,j+1)+ &
947 74480 : cofum_solv(i,j,4)*phisolv(i-1,j+1)+ &
948 : cofum_solv(i,j,5)*phisolv(i-1,j)+ &
949 74480 : cofum_solv(i,j,6)*phisolv(i-1,j-1)+ &
950 : cofum_solv(i,j,7)*phisolv(i,j-1)+ &
951 : cofum_solv(i,j,8)*phisolv(i+1,j-1)+ &
952 223440 : cofum_solv(i,j,9)*phisolv(i,j))
953 :
954 74480 : ressolv(i,j) = c0(i,j,10)-ressolv(i,j)
955 75411 : l2norm = l2norm + ressolv(i,j)*ressolv(i,j)
956 : enddo
957 : enddo
958 : ! write(iulog,*) 'L2norm (global root task) = ',l2norm
959 :
960 19 : end subroutine solver
961 : !--------------------------------------------------------------------
962 : end module edyn_solve
|
