31 real,
parameter::
r3 = 1./3.
35 real,
parameter::
a1 = 0.5625
36 real,
parameter::
a2 = -0.0625
40 real,
parameter::
b1 = 7./12.
41 real,
parameter::
b2 = -1./12.
48 #ifndef USE_OLD_ALGORITHM 49 subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
50 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
51 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
52 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
54 logical,
optional,
intent(IN):: replace
56 real,
parameter:: c1 = 2./3.
57 real,
parameter:: c2 = -1./6.
62 real qx(is:ie+1,js-ng:je+ng)
63 real qy(is-ng:ie+ng,js:je+1)
64 real qxx(is-ng:ie+ng,js-ng:je+ng)
65 real qyy(is-ng:ie+ng,js-ng:je+ng)
68 real:: q1(is-1:ie+1), q2(js-1:je+1)
69 integer:: i, j, is1, js1, is2, js2, ie1, je1
71 real,
pointer,
dimension(:,:,:) :: grid, agrid
72 real,
pointer,
dimension(:,:) :: dxa, dya
73 real(kind=R_GRID),
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
75 edge_w => gridstruct%edge_w
76 edge_e => gridstruct%edge_e
77 edge_s => gridstruct%edge_s
78 edge_n => gridstruct%edge_n
80 grid => gridstruct%grid
81 agrid => gridstruct%agrid
85 if (gridstruct%grid_type < 3)
then 97 if (gridstruct%nested)
then 101 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
108 if ( gridstruct%sw_corner )
then 109 p0(1:2) = grid(1,1,1:2)
110 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
111 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
112 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 115 if ( gridstruct%se_corner )
then 116 p0(1:2) = grid(npx,1,1:2)
117 qout(npx,1) = (
extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
118 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
119 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 121 if ( gridstruct%ne_corner )
then 122 p0(1:2) = grid(npx,npy,1:2)
123 qout(npx,npy) = (
extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
124 extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
125 extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*
r3 127 if ( gridstruct%nw_corner )
then 128 p0(1:2) = grid(1,npy,1:2)
129 qout(1,npy) = (
extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
130 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
131 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 137 do j=
max(1,js-2),
min(npy-1,je+2)
138 do i=
max(3,is),
min(npx-2,ie+1)
139 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
146 q2(j) = (qin(0,j)*dxa(1,j) + qin(1,j)*dxa(0,j))/(dxa(0,j) + dxa(1,j))
149 qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
152 do j=
max(1,js-2),
min(npy-1,je+2)
153 g_in = dxa(2,j) / dxa(1,j)
154 g_ou = dxa(-1,j) / dxa(0,j)
155 qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
156 ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
157 qx(2,j) = ( 3.*(g_in*qin(1,j)+qin(2,j))-(g_in*qx(1,j)+qx(3,j)) ) / (2.+2.*g_in)
162 if ( (ie+1)==npx )
then 164 q2(j) = (qin(npx-1,j)*dxa(npx,j) + qin(npx,j)*dxa(npx-1,j))/(dxa(npx-1,j) + dxa(npx,j))
167 qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
170 do j=
max(1,js-2),
min(npy-1,je+2)
171 g_in = dxa(npx-2,j) / dxa(npx-1,j)
172 g_ou = dxa(npx+1,j) / dxa(npx,j)
173 qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
174 ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
175 qx(npx-1,j) = (3.*(qin(npx-2,j)+g_in*qin(npx-1,j)) - (g_in*qx(npx,j)+qx(npx-2,j)))/(2.+2.*g_in)
184 if (gridstruct%nested)
then 189 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1) + qin(i,j))
195 do j=
max(3,js),
min(npy-2,je+1)
196 do i=
max(1,is-2),
min(npx-1,ie+2)
197 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1) + qin(i,j))
204 q1(i) = (qin(i,0)*dya(i,1) + qin(i,1)*dya(i,0))/(dya(i,0) + dya(i,1))
207 qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
210 do i=
max(1,is-2),
min(npx-1,ie+2)
211 g_in = dya(i,2) / dya(i,1)
212 g_ou = dya(i,-1) / dya(i,0)
213 qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
214 ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
215 qy(i,2) = (3.*(g_in*qin(i,1)+qin(i,2)) - (g_in*qy(i,1)+qy(i,3)))/(2.+2.*g_in)
220 if ( (je+1)==npy )
then 222 q1(i) = (qin(i,npy-1)*dya(i,npy) + qin(i,npy)*dya(i,npy-1))/(dya(i,npy-1)+dya(i,npy))
225 qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
228 do i=
max(1,is-2),
min(npx-1,ie+2)
229 g_in = dya(i,npy-2) / dya(i,npy-1)
230 g_ou = dya(i,npy+1) / dya(i,npy)
231 qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
232 ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
233 qy(i,npy-1) = (3.*(qin(i,npy-2)+g_in*qin(i,npy-1)) - (g_in*qy(i,npy)+qy(i,npy-2)))/(2.+2.*g_in)
240 if (gridstruct%nested)
then 244 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
250 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
254 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
262 do j=
max(3,js),
min(npy-2,je+1)
263 do i=
max(2,is),
min(npx-1,ie+1)
264 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
269 do i=
max(2,is),
min(npx-1,ie+1)
270 qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
273 if ( (je+1)==npy )
then 274 do i=
max(2,is),
min(npx-1,ie+1)
275 qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
280 do j=
max(2,js),
min(npy-1,je+1)
281 do i=
max(3,is),
min(npx-2,ie+1)
282 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
284 if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
285 if((ie+1)==npx) qyy(npx-1,j) = c1*(qy(npx-2,j)+qy(npx-1,j))+c2*(qout(npx,j)+qyy(npx-2,j))
287 do i=
max(2,is),
min(npx-1,ie+1)
288 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
301 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
307 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
313 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
314 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
319 if (
present(replace) )
then 334 subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
335 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
336 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
337 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
339 logical,
optional,
intent(IN):: replace
341 real,
parameter:: c1 = 2./3.
342 real,
parameter:: c2 = -1./6.
349 real,
parameter:: d1 = 0.375
350 real,
parameter:: d2 = -1./24.
352 real qx(is:ie+1,js-ng:je+ng)
353 real qy(is-ng:ie+ng,js:je+1)
354 real qxx(is-ng:ie+ng,js-ng:je+ng)
355 real qyy(is-ng:ie+ng,js-ng:je+ng)
359 real q1(npx), q2(npy)
360 integer:: i, j, is1, js1, is2, js2, ie1, je1
362 real,
pointer,
dimension(:,:,:) :: grid, agrid
363 real,
pointer,
dimension(:,:) :: dxa, dya
364 real,
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
366 edge_w => gridstruct%edge_w
367 edge_e => gridstruct%edge_e
368 edge_s => gridstruct%edge_s
369 edge_n => gridstruct%edge_n
372 grid => gridstruct%grid
373 agrid => gridstruct%agrid
374 dxa => gridstruct%dxa
375 dya => gridstruct%dya
377 if (gridstruct%grid_type < 3)
then 384 ie1 =
min(npx-1,ie+1)
385 je1 =
min(npy-1,je+1)
389 if ( gridstruct%sw_corner ) qout(1, 1) =
r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
390 if ( gridstruct%se_corner ) qout(npx, 1) =
r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
391 if ( gridstruct%ne_corner ) qout(npx,npy) =
r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
392 if ( gridstruct%nw_corner ) qout(1, npy) =
r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
397 if ( gridstruct%sw_corner )
then 398 qout(1,1) = d1*(qin(1, 0) + qin( 0,1) + qin(1,1)) + &
399 d2*(qin(2,-1) + qin(-1,2) + qin(2,2))
401 if ( gridstruct%se_corner )
then 402 qout(npx,1) = d1*(qin(npx-1, 0) + qin(npx-1,1) + qin(npx, 1)) + &
403 d2*(qin(npx-2,-1) + qin(npx-2,2) + qin(npx+1,2))
405 if ( gridstruct%ne_corner )
then 406 qout(npx,npy) = d1*(qin(npx-1,npy-1) + qin(npx, npy-1) + qin(npx-1,npy)) + &
407 d2*(qin(npx-2,npy-2) + qin(npx+1,npy-2) + qin(npx-2,npy+1))
409 if ( gridstruct%nw_corner )
then 410 qout(1,npy) = d1*(qin( 0,npy-1) + qin(1,npy-1) + qin(1,npy)) + &
411 d2*(qin(-1,npy-2) + qin(2,npy-2) + qin(2,npy+1))
415 if ( gridstruct%sw_corner )
then 416 p0(1:2) = grid(1,1,1:2)
417 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
418 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
419 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 422 if ( gridstruct%se_corner )
then 423 p0(1:2) = grid(npx,1,1:2)
424 qout(npx,1) = (
extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
425 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
426 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 428 if ( gridstruct%ne_corner )
then 429 p0(1:2) = grid(npx,npy,1:2)
430 qout(npx,npy) = (
extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
431 extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
432 extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*
r3 434 if ( gridstruct%nw_corner )
then 435 p0(1:2) = grid(1,npy,1:2)
436 qout(1,npy) = (
extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
437 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
438 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 446 if (gridstruct%nested)
then 450 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
457 do j=
max(1,js-2),
min(npy-1,je+2)
458 do i=
max(3,is),
min(npx-2,ie+1)
459 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
466 do j=
max(1,js-2),
min(npy-1,je+2)
467 gratio = dxa(2,j) / dxa(1,j)
469 qx(1,j) = 0.5*((2.+gratio)*(qin(0,j)+qin(1,j))-(qin(-1,j)+qin(2,j))) / (1.+gratio)
472 g_ou = dxa(-1,j) / dxa(0,j)
473 qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
474 ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
476 qx(2,j) = ( 3.*(gratio*qin(1,j)+qin(2,j))-(gratio*qx(1,j)+qx(3,j)) ) / (2.+2.*gratio)
479 do j=
max(3,js),
min(npy-2,je+1)
480 qout(1,j) =
a2*(qx(1,j-2)+qx(1,j+1)) +
a1*(qx(1,j-1)+qx(1,j))
483 if( js==1 ) qout(1, 2) = c1*(qx(1,1)+qx(1,2)) + c2*(qout(1,1)+qout(1,3))
484 if((je+1)==npy) qout(1,npy-1) = c1*(qx(1,npy-2)+qx(1,npy-1)) + c2*(qout(1,npy-2)+qout(1,npy))
488 if ( (ie+1)==npx )
then 490 do j=
max(1,js-2),
min(npy-1,je+2)
491 gratio = dxa(npx-2,j) / dxa(npx-1,j)
493 qx(npx,j) = 0.5*((2.+gratio)*(qin(npx-1,j)+qin(npx,j)) &
494 - (qin(npx-2,j)+qin(npx+1,j))) / (1.+gratio )
497 g_ou = dxa(npx+1,j) / dxa(npx,j)
498 qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
499 ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
501 qx(npx-1,j) = (3.*(qin(npx-2,j)+gratio*qin(npx-1,j)) - (gratio*qx(npx,j)+qx(npx-2,j)))/(2.+2.*gratio)
504 do j=
max(3,js),
min(npy-2,je+1)
505 qout(npx,j) =
a2*(qx(npx,j-2)+qx(npx,j+1)) +
a1*(qx(npx,j-1)+qx(npx,j))
508 if(js==1) qout(npx,2) = c1*(qx(npx,1)+qx(npx,2))+c2*(qout(npx,1)+qout(npx,3))
509 if((je+1)==npy) qout(npx,npy-1) = &
510 c1*(qx(npx,npy-2)+qx(npx,npy-1))+c2*(qout(npx,npy-2)+qout(npx,npy))
517 if (gridstruct%nested)
then 521 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1)+qin(i,j))
527 do j=
max(3,js),
min(npy-2,je+1)
528 do i=
max(1,is-2),
min(npx-1,ie+2)
529 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1)+qin(i,j))
536 do i=
max(1,is-2),
min(npx-1,ie+2)
537 gratio = dya(i,2) / dya(i,1)
539 qy(i,1) = 0.5*((2.+gratio)*(qin(i,0)+qin(i,1)) &
540 - (qin(i,-1)+qin(i,2))) / (1.+gratio )
543 g_ou = dya(i,-1) / dya(i,0)
544 qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
545 ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
547 qy(i,2) = (3.*(gratio*qin(i,1)+qin(i,2)) - (gratio*qy(i,1)+qy(i,3)))/(2.+2.*gratio)
550 do i=
max(3,is),
min(npx-2,ie+1)
551 qout(i,1) =
a2*(qy(i-2,1)+qy(i+1,1)) +
a1*(qy(i-1,1)+qy(i,1))
554 if( is==1 ) qout(2,1) = c1*(qy(1,1)+qy(2,1))+c2*(qout(1,1)+qout(3,1))
555 if((ie+1)==npx) qout(npx-1,1) = c1*(qy(npx-2,1)+qy(npx-1,1))+c2*(qout(npx-2,1)+qout(npx,1))
560 if ( (je+1)==npy )
then 561 do i=
max(1,is-2),
min(npx-1,ie+2)
562 gratio = dya(i,npy-2) / dya(i,npy-1)
564 qy(i,npy ) = 0.5*((2.+gratio)*(qin(i,npy-1)+qin(i,npy)) &
565 - (qin(i,npy-2)+qin(i,npy+1))) / (1.+gratio)
568 g_ou = dya(i,npy+1) / dya(i,npy)
569 qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
570 ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
572 qy(i,npy-1) = (3.*(qin(i,npy-2)+gratio*qin(i,npy-1)) - (gratio*qy(i,npy)+qy(i,npy-2)))/(2.+2.*gratio)
575 do i=
max(3,is),
min(npx-2,ie+1)
576 qout(i,npy) =
a2*(qy(i-2,npy)+qy(i+1,npy)) +
a1*(qy(i-1,npy)+qy(i,npy))
579 if( is==1 ) qout(2,npy) = c1*(qy(1,npy)+qy(2,npy))+c2*(qout(1,npy)+qout(3,npy))
580 if((ie+1)==npx) qout(npx-1,npy) = c1*(qy(npx-2,npy)+qy(npx-1,npy))+c2*(qout(npx-2,npy)+qout(npx,npy))
585 if (gridstruct%nested)
then 589 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
595 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
599 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
606 do j=
max(3,js),
min(npy-2,je+1)
607 do i=
max(2,is),
min(npx-1,ie+1)
608 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
613 do i=
max(2,is),
min(npx-1,ie+1)
614 qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
617 if ( (je+1)==npy )
then 618 do i=
max(2,is),
min(npx-1,ie+1)
619 qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
624 do j=
max(2,js),
min(npy-1,je+1)
625 do i=
max(3,is),
min(npx-2,ie+1)
626 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
628 if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
629 if((ie+1)==npx) qyy(npx-1,j) = c1*(qy(npx-2,j)+qy(npx-1,j))+c2*(qout(npx,j)+qyy(npx-2,j))
631 do i=
max(2,is),
min(npx-1,ie+1)
632 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
645 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
651 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
657 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
658 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
663 if (
present(replace) )
then 677 subroutine a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
678 integer,
intent(IN ) :: npx, npy, is, ie, js, je, ng
679 real ,
intent(INOUT) :: qin(is-ng:ie+ng,js-ng:je+ng)
680 real ,
intent( OUT) :: qout(is-ng:ie+ng,js-ng:je+ng)
682 logical,
optional,
intent(IN) :: replace
684 real q1(npx), q2(npy)
686 integer :: is1, js1, is2, js2, ie1, je1
688 real,
pointer,
dimension(:,:,:) :: grid, agrid
689 real,
pointer,
dimension(:,:) :: dxa, dya
691 real(kind=R_GRID),
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
693 edge_w => gridstruct%edge_w
694 edge_e => gridstruct%edge_e
695 edge_s => gridstruct%edge_s
696 edge_n => gridstruct%edge_n
698 grid => gridstruct%grid
699 agrid => gridstruct%agrid
700 dxa => gridstruct%dxa
701 dya => gridstruct%dya
703 if (gridstruct%grid_type < 3)
then 705 if (gridstruct%nested)
then 709 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
720 ie1 =
min(npx-1,ie+1)
721 je1 =
min(npy-1,je+1)
725 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
730 if ( gridstruct%sw_corner ) qout(1, 1) =
r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
731 if ( gridstruct%se_corner ) qout(npx, 1) =
r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
732 if ( gridstruct%ne_corner ) qout(npx,npy) =
r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
733 if ( gridstruct%nw_corner ) qout(1, npy) =
r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
738 q2(j) = 0.5*(qin(0,j) + qin(1,j))
741 qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
746 if ( (ie+1)==npx )
then 748 q2(j) = 0.5*(qin(npx-1,j) + qin(npx,j))
751 qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
758 q1(i) = 0.5*(qin(i,0) + qin(i,1))
761 qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
766 if ( (je+1)==npy )
then 768 q1(i) = 0.5*(qin(i,npy-1) + qin(i,npy))
771 qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
781 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
788 if (
present(replace) )
then 801 real,
intent(in ),
dimension(2):: p0, p1, p2
802 real,
intent(in ):: q1, q2
805 x1 = great_circle_dist(
real(p1,kind=R_GRID),
real(p0,kind=R_GRID) )
806 x2 = great_circle_dist(
real(p2,kind=R_GRID),
real(p0,kind=R_GRID) )
813 subroutine a2b_ord4(qin, qout, grid, agrid, npx, npy, is, ie, js, je, ng, replace)
815 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
816 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
817 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
818 real,
intent(in) :: grid(is-ng:ie+ng+1,js-ng:je+ng+1,2)
819 real,
intent(in) :: agrid(is-ng:ie+ng,js-ng:je+ng,2)
820 logical,
optional,
intent(IN):: replace
821 real qx(is:ie+1,js-ng:je+ng)
822 real qy(is-ng:ie+ng,js:je+1)
826 real,
pointer,
dimension(:,:,:) :: grid, agrid
827 real,
pointer,
dimension(:,:) :: dxa, dya
829 real,
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
831 edge_w => gridstruct%edge_w
832 edge_e => gridstruct%edge_e
833 edge_s => gridstruct%edge_s
834 edge_n => gridstruct%edge_n
836 grid => gridstruct%grid
837 agrid => gridstruct%agrid
838 dxa => gridstruct%dxa
839 dya => gridstruct%dya
842 if (gridstruct%grid_type < 3)
then 852 if ( i==1 .and. j==1 )
goto 123
853 if ( i==2 .and. j==1 )
then 854 qin(0,-1) = qin(-1,2)
855 qin(0, 0) = qin(-1,1)
857 if ( i==1 .and. j==2 )
then 858 qin(-1,0) = qin(2,-1)
859 qin( 0,0) = qin(1,-1)
861 if ( i==2 .and. j==2 )
then 865 if ( i==npx .and. j==1 )
goto 123
866 if ( i==npx-1 .and. j==1 )
then 867 qin(npx,-1) = qin(npx+1,2)
868 qin(npx, 0) = qin(npx+1,1)
870 if ( i==npx-1 .and. j==2 )
then 871 qin(npx,0) = qin(npx-4,4)
873 if ( i==npx .and. j==2 )
then 874 qin(npx+1,0) = qin(npx-2,-1)
875 qin(npx, 0) = qin(npx-1,-1)
878 if ( i==npx .and. j==npy )
goto 123
879 if ( i==npx-1 .and. j==npy-1 )
then 880 qin(npx,npy) = qin(npx-4,npy-4)
882 if ( i==npx .and. j==npy-1 )
then 883 qin(npx+1,npy) = qin(npx-2,npy+1)
884 qin(npx, npy) = qin(npx-1,npy+1)
886 if ( i==npx-1 .and. j==npy )
then 887 qin(npx,npy+1) = qin(npx+1,npy-2)
888 qin(npx,npy ) = qin(npx+1,npy-1)
891 if ( i==1 .and. j==npy )
goto 123
892 if ( i==1 .and. j==npy-1 )
then 893 qin(-1,npy) = qin(2,npy+1)
894 qin( 0,npy) = qin(1,npy+1)
896 if ( i==2 .and. j==npy-1 )
then 897 qin(0,npy) = qin(4,npy-4)
899 if ( i==2 .and. j==npy )
then 900 qin(0,npy+1) = qin(-1,npy-2)
901 qin(0,npy ) = qin(-1,npy-1)
904 qout(i,j) = van2(1, i,j)*qin(i-2,j-2) + van2(2, i,j)*qin(i-1,j-2) + &
905 van2(3, i,j)*qin(i ,j-2) + van2(4, i,j)*qin(i+1,j-2) + &
906 van2(5, i,j)*qin(i-2,j-1) + van2(6, i,j)*qin(i-1,j-1) + &
907 van2(7, i,j)*qin(i ,j-1) + van2(8, i,j)*qin(i+1,j-1) + &
908 van2(9, i,j)*qin(i-2,j ) + van2(10,i,j)*qin(i-1,j ) + &
909 van2(11,i,j)*qin(i ,j ) + van2(12,i,j)*qin(i+1,j ) + &
910 van2(13,i,j)*qin(i-2,j+1) + van2(14,i,j)*qin(i-1,j+1) + &
911 van2(15,i,j)*qin(i ,j+1) + van2(16,i,j)*qin(i+1,j+1)
917 if ( gridstruct%sw_corner )
then 918 p0(1:2) = grid(1,1,1:2)
919 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
920 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
921 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 924 if ( gridstruct%se_corner )
then 925 p0(1:2) = grid(npx,1,1:2)
926 qout(npx,1) = (
extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
927 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
928 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 930 if ( gridstruct%ne_corner )
then 931 p0(1:2) = grid(npx,npy,1:2)
932 qout(npx,npy) = (
extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
933 extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
934 extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*
r3 936 if ( gridstruct%nw_corner )
then 937 p0(1:2) = grid(1,npy,1:2)
938 qout(1,npy) = (
extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
939 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
940 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 951 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
957 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
963 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
964 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
970 if (
present(replace) )
then
void a2b_ord2(int nx, int ny, const double *qin, const double *edge_w, const double *edge_e, const double *edge_s, const double *edge_n, double *qout, int on_west_edge, int on_east_edge, int on_south_edge, int on_north_edge)
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
integer, parameter, public r_grid
real function, public great_circle_dist(q1, q2, radius)
real function extrap_corner(p0, p1, p2, q1, q2)