51 real,
parameter::
r3 = 1./3.
55 real,
parameter::
a1 = 0.5625
56 real,
parameter::
a2 = -0.0625
60 real,
parameter::
b1 = 7./12.
61 real,
parameter::
b2 = -1./12.
68 #ifndef USE_OLD_ALGORITHM 69 subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
70 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
71 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
72 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
74 logical,
optional,
intent(IN):: replace
76 real,
parameter:: c1 = 2./3.
77 real,
parameter:: c2 = -1./6.
82 real qx(is:ie+1,js-ng:je+ng)
83 real qy(is-ng:ie+ng,js:je+1)
84 real qxx(is-ng:ie+ng,js-ng:je+ng)
85 real qyy(is-ng:ie+ng,js-ng:je+ng)
88 real:: q1(is-1:ie+1), q2(js-1:je+1)
89 integer:: i, j, is1, js1, is2, js2, ie1, je1
91 real,
pointer,
dimension(:,:,:) :: grid, agrid
92 real,
pointer,
dimension(:,:) :: dxa, dya
93 real(kind=R_GRID),
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
95 edge_w => gridstruct%edge_w
96 edge_e => gridstruct%edge_e
97 edge_s => gridstruct%edge_s
98 edge_n => gridstruct%edge_n
100 grid => gridstruct%grid
101 agrid => gridstruct%agrid
102 dxa => gridstruct%dxa
103 dya => gridstruct%dya
105 if (gridstruct%grid_type < 3)
then 112 ie1 = min(npx-1,ie+1)
113 je1 = min(npy-1,je+1)
117 if (gridstruct%bounded_domain)
then 121 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
128 if ( gridstruct%sw_corner )
then 129 p0(1:2) = grid(1,1,1:2)
130 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
131 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
132 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 135 if ( gridstruct%se_corner )
then 136 p0(1:2) = grid(npx,1,1:2)
137 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)) + &
138 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
139 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 141 if ( gridstruct%ne_corner )
then 142 p0(1:2) = grid(npx,npy,1:2)
143 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)) + &
144 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)) + &
145 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 147 if ( gridstruct%nw_corner )
then 148 p0(1:2) = grid(1,npy,1:2)
149 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)) + &
150 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
151 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 157 do j=max(1,js-2),min(npy-1,je+2)
158 do i=max(3,is), min(npx-2,ie+1)
159 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
166 q2(j) = (qin(0,j)*dxa(1,j) + qin(1,j)*dxa(0,j))/(dxa(0,j) + dxa(1,j))
169 qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
172 do j=max(1,js-2),min(npy-1,je+2)
173 g_in = dxa(2,j) / dxa(1,j)
174 g_ou = dxa(-1,j) / dxa(0,j)
175 qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
176 ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
177 qx(2,j) = ( 3.*(g_in*qin(1,j)+qin(2,j))-(g_in*qx(1,j)+qx(3,j)) ) / (2.+2.*g_in)
182 if ( (ie+1)==npx )
then 184 q2(j) = (qin(npx-1,j)*dxa(npx,j) + qin(npx,j)*dxa(npx-1,j))/(dxa(npx-1,j) + dxa(npx,j))
187 qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
190 do j=max(1,js-2),min(npy-1,je+2)
191 g_in = dxa(npx-2,j) / dxa(npx-1,j)
192 g_ou = dxa(npx+1,j) / dxa(npx,j)
193 qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
194 ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
195 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)
204 if (gridstruct%bounded_domain)
then 209 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1) + qin(i,j))
215 do j=max(3,js),min(npy-2,je+1)
216 do i=max(1,is-2), min(npx-1,ie+2)
217 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1) + qin(i,j))
224 q1(i) = (qin(i,0)*dya(i,1) + qin(i,1)*dya(i,0))/(dya(i,0) + dya(i,1))
227 qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
230 do i=max(1,is-2),min(npx-1,ie+2)
231 g_in = dya(i,2) / dya(i,1)
232 g_ou = dya(i,-1) / dya(i,0)
233 qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
234 ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
235 qy(i,2) = (3.*(g_in*qin(i,1)+qin(i,2)) - (g_in*qy(i,1)+qy(i,3)))/(2.+2.*g_in)
240 if ( (je+1)==npy )
then 242 q1(i) = (qin(i,npy-1)*dya(i,npy) + qin(i,npy)*dya(i,npy-1))/(dya(i,npy-1)+dya(i,npy))
245 qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
248 do i=max(1,is-2),min(npx-1,ie+2)
249 g_in = dya(i,npy-2) / dya(i,npy-1)
250 g_ou = dya(i,npy+1) / dya(i,npy)
251 qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
252 ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
253 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)
260 if (gridstruct%bounded_domain)
then 264 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
270 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
274 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
282 do j=max(3,js),min(npy-2,je+1)
283 do i=max(2,is),min(npx-1,ie+1)
284 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
289 do i=max(2,is),min(npx-1,ie+1)
290 qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
293 if ( (je+1)==npy )
then 294 do i=max(2,is),min(npx-1,ie+1)
295 qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
300 do j=max(2,js),min(npy-1,je+1)
301 do i=max(3,is),min(npx-2,ie+1)
302 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
304 if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
305 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))
307 do i=max(2,is),min(npx-1,ie+1)
308 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
321 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
327 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
333 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
334 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
339 if (
present(replace) )
then 354 subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
355 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
356 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
357 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
359 logical,
optional,
intent(IN):: replace
361 real,
parameter:: c1 = 2./3.
362 real,
parameter:: c2 = -1./6.
369 real,
parameter:: d1 = 0.375
370 real,
parameter:: d2 = -1./24.
372 real qx(is:ie+1,js-ng:je+ng)
373 real qy(is-ng:ie+ng,js:je+1)
374 real qxx(is-ng:ie+ng,js-ng:je+ng)
375 real qyy(is-ng:ie+ng,js-ng:je+ng)
379 real q1(npx), q2(npy)
380 integer:: i, j, is1, js1, is2, js2, ie1, je1
382 real,
pointer,
dimension(:,:,:) :: grid, agrid
383 real,
pointer,
dimension(:,:) :: dxa, dya
384 real,
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
386 edge_w => gridstruct%edge_w
387 edge_e => gridstruct%edge_e
388 edge_s => gridstruct%edge_s
389 edge_n => gridstruct%edge_n
392 grid => gridstruct%grid
393 agrid => gridstruct%agrid
394 dxa => gridstruct%dxa
395 dya => gridstruct%dya
397 if (gridstruct%grid_type < 3)
then 404 ie1 = min(npx-1,ie+1)
405 je1 = min(npy-1,je+1)
409 if ( gridstruct%sw_corner ) qout(1, 1) =
r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
410 if ( gridstruct%se_corner ) qout(npx, 1) =
r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
411 if ( gridstruct%ne_corner ) qout(npx,npy) =
r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
412 if ( gridstruct%nw_corner ) qout(1, npy) =
r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
417 if ( gridstruct%sw_corner )
then 418 qout(1,1) = d1*(qin(1, 0) + qin( 0,1) + qin(1,1)) + &
419 d2*(qin(2,-1) + qin(-1,2) + qin(2,2))
421 if ( gridstruct%se_corner )
then 422 qout(npx,1) = d1*(qin(npx-1, 0) + qin(npx-1,1) + qin(npx, 1)) + &
423 d2*(qin(npx-2,-1) + qin(npx-2,2) + qin(npx+1,2))
425 if ( gridstruct%ne_corner )
then 426 qout(npx,npy) = d1*(qin(npx-1,npy-1) + qin(npx, npy-1) + qin(npx-1,npy)) + &
427 d2*(qin(npx-2,npy-2) + qin(npx+1,npy-2) + qin(npx-2,npy+1))
429 if ( gridstruct%nw_corner )
then 430 qout(1,npy) = d1*(qin( 0,npy-1) + qin(1,npy-1) + qin(1,npy)) + &
431 d2*(qin(-1,npy-2) + qin(2,npy-2) + qin(2,npy+1))
435 if ( gridstruct%sw_corner )
then 436 p0(1:2) = grid(1,1,1:2)
437 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
438 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
439 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 442 if ( gridstruct%se_corner )
then 443 p0(1:2) = grid(npx,1,1:2)
444 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)) + &
445 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
446 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 448 if ( gridstruct%ne_corner )
then 449 p0(1:2) = grid(npx,npy,1:2)
450 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)) + &
451 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)) + &
452 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 454 if ( gridstruct%nw_corner )
then 455 p0(1:2) = grid(1,npy,1:2)
456 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)) + &
457 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
458 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 466 if (gridstruct%bounded_domain)
then 470 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
477 do j=max(1,js-2),min(npy-1,je+2)
478 do i=max(3,is), min(npx-2,ie+1)
479 qx(i,j) =
b2*(qin(i-2,j)+qin(i+1,j)) +
b1*(qin(i-1,j)+qin(i,j))
486 do j=max(1,js-2),min(npy-1,je+2)
487 gratio = dxa(2,j) / dxa(1,j)
489 qx(1,j) = 0.5*((2.+gratio)*(qin(0,j)+qin(1,j))-(qin(-1,j)+qin(2,j))) / (1.+gratio)
492 g_ou = dxa(-1,j) / dxa(0,j)
493 qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
494 ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
496 qx(2,j) = ( 3.*(gratio*qin(1,j)+qin(2,j))-(gratio*qx(1,j)+qx(3,j)) ) / (2.+2.*gratio)
499 do j=max(3,js),min(npy-2,je+1)
500 qout(1,j) =
a2*(qx(1,j-2)+qx(1,j+1)) +
a1*(qx(1,j-1)+qx(1,j))
503 if( js==1 ) qout(1, 2) = c1*(qx(1,1)+qx(1,2)) + c2*(qout(1,1)+qout(1,3))
504 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))
508 if ( (ie+1)==npx )
then 510 do j=max(1,js-2),min(npy-1,je+2)
511 gratio = dxa(npx-2,j) / dxa(npx-1,j)
513 qx(npx,j) = 0.5*((2.+gratio)*(qin(npx-1,j)+qin(npx,j)) &
514 - (qin(npx-2,j)+qin(npx+1,j))) / (1.+gratio )
517 g_ou = dxa(npx+1,j) / dxa(npx,j)
518 qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
519 ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
521 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)
524 do j=max(3,js),min(npy-2,je+1)
525 qout(npx,j) =
a2*(qx(npx,j-2)+qx(npx,j+1)) +
a1*(qx(npx,j-1)+qx(npx,j))
528 if(js==1) qout(npx,2) = c1*(qx(npx,1)+qx(npx,2))+c2*(qout(npx,1)+qout(npx,3))
529 if((je+1)==npy) qout(npx,npy-1) = &
530 c1*(qx(npx,npy-2)+qx(npx,npy-1))+c2*(qout(npx,npy-2)+qout(npx,npy))
537 if (gridstruct%bounded_domain)
then 541 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1)+qin(i,j))
547 do j=max(3,js),min(npy-2,je+1)
548 do i=max(1,is-2), min(npx-1,ie+2)
549 qy(i,j) =
b2*(qin(i,j-2)+qin(i,j+1)) +
b1*(qin(i,j-1)+qin(i,j))
556 do i=max(1,is-2),min(npx-1,ie+2)
557 gratio = dya(i,2) / dya(i,1)
559 qy(i,1) = 0.5*((2.+gratio)*(qin(i,0)+qin(i,1)) &
560 - (qin(i,-1)+qin(i,2))) / (1.+gratio )
563 g_ou = dya(i,-1) / dya(i,0)
564 qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
565 ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
567 qy(i,2) = (3.*(gratio*qin(i,1)+qin(i,2)) - (gratio*qy(i,1)+qy(i,3)))/(2.+2.*gratio)
570 do i=max(3,is),min(npx-2,ie+1)
571 qout(i,1) =
a2*(qy(i-2,1)+qy(i+1,1)) +
a1*(qy(i-1,1)+qy(i,1))
574 if( is==1 ) qout(2,1) = c1*(qy(1,1)+qy(2,1))+c2*(qout(1,1)+qout(3,1))
575 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))
580 if ( (je+1)==npy )
then 581 do i=max(1,is-2),min(npx-1,ie+2)
582 gratio = dya(i,npy-2) / dya(i,npy-1)
584 qy(i,npy ) = 0.5*((2.+gratio)*(qin(i,npy-1)+qin(i,npy)) &
585 - (qin(i,npy-2)+qin(i,npy+1))) / (1.+gratio)
588 g_ou = dya(i,npy+1) / dya(i,npy)
589 qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
590 ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
592 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)
595 do i=max(3,is),min(npx-2,ie+1)
596 qout(i,npy) =
a2*(qy(i-2,npy)+qy(i+1,npy)) +
a1*(qy(i-1,npy)+qy(i,npy))
599 if( is==1 ) qout(2,npy) = c1*(qy(1,npy)+qy(2,npy))+c2*(qout(1,npy)+qout(3,npy))
600 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))
605 if (gridstruct%bounded_domain)
then 609 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
615 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
619 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
626 do j=max(3,js),min(npy-2,je+1)
627 do i=max(2,is),min(npx-1,ie+1)
628 qxx(i,j) =
a2*(qx(i,j-2)+qx(i,j+1)) +
a1*(qx(i,j-1)+qx(i,j))
633 do i=max(2,is),min(npx-1,ie+1)
634 qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
637 if ( (je+1)==npy )
then 638 do i=max(2,is),min(npx-1,ie+1)
639 qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
644 do j=max(2,js),min(npy-1,je+1)
645 do i=max(3,is),min(npx-2,ie+1)
646 qyy(i,j) =
a2*(qy(i-2,j)+qy(i+1,j)) +
a1*(qy(i-1,j)+qy(i,j))
648 if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
649 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))
651 do i=max(2,is),min(npx-1,ie+1)
652 qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j))
665 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
671 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
677 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
678 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
683 if (
present(replace) )
then 697 subroutine a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
698 integer,
intent(IN ) :: npx, npy, is, ie, js, je, ng
699 real ,
intent(INOUT) :: qin(is-ng:ie+ng,js-ng:je+ng)
700 real ,
intent( OUT) :: qout(is-ng:ie+ng,js-ng:je+ng)
702 logical,
optional,
intent(IN) :: replace
704 real q1(npx), q2(npy)
706 integer :: is1, js1, is2, js2, ie1, je1
708 real,
pointer,
dimension(:,:,:) :: grid, agrid
709 real,
pointer,
dimension(:,:) :: dxa, dya
711 real(kind=R_GRID),
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
713 edge_w => gridstruct%edge_w
714 edge_e => gridstruct%edge_e
715 edge_s => gridstruct%edge_s
716 edge_n => gridstruct%edge_n
718 grid => gridstruct%grid
719 agrid => gridstruct%agrid
720 dxa => gridstruct%dxa
721 dya => gridstruct%dya
723 if (gridstruct%grid_type < 3)
then 725 if (gridstruct%bounded_domain)
then 729 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
740 ie1 = min(npx-1,ie+1)
741 je1 = min(npy-1,je+1)
745 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
750 if ( gridstruct%sw_corner ) qout(1, 1) =
r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
751 if ( gridstruct%se_corner ) qout(npx, 1) =
r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
752 if ( gridstruct%ne_corner ) qout(npx,npy) =
r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
753 if ( gridstruct%nw_corner ) qout(1, npy) =
r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
758 q2(j) = 0.5*(qin(0,j) + qin(1,j))
761 qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
766 if ( (ie+1)==npx )
then 768 q2(j) = 0.5*(qin(npx-1,j) + qin(npx,j))
771 qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
778 q1(i) = 0.5*(qin(i,0) + qin(i,1))
781 qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
786 if ( (je+1)==npy )
then 788 q1(i) = 0.5*(qin(i,npy-1) + qin(i,npy))
791 qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
801 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
808 if (
present(replace) )
then 821 real,
intent(in ),
dimension(2):: p0, p1, p2
822 real,
intent(in ):: q1, q2
825 x1 = great_circle_dist(
real(p1,kind=R_GRID),
real(p0,kind=R_GRID) )
826 x2 = great_circle_dist(
real(p2,kind=R_GRID),
real(p0,kind=R_GRID) )
833 subroutine a2b_ord4(qin, qout, grid, agrid, npx, npy, is, ie, js, je, ng, replace)
835 integer,
intent(IN):: npx, npy, is, ie, js, je, ng
836 real,
intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng)
837 real,
intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng)
838 real,
intent(in) :: grid(is-ng:ie+ng+1,js-ng:je+ng+1,2)
839 real,
intent(in) :: agrid(is-ng:ie+ng,js-ng:je+ng,2)
840 logical,
optional,
intent(IN):: replace
841 real qx(is:ie+1,js-ng:je+ng)
842 real qy(is-ng:ie+ng,js:je+1)
846 real,
pointer,
dimension(:,:,:) :: grid, agrid
847 real,
pointer,
dimension(:,:) :: dxa, dya
849 real,
pointer,
dimension(:) :: edge_w, edge_e, edge_s, edge_n
851 edge_w => gridstruct%edge_w
852 edge_e => gridstruct%edge_e
853 edge_s => gridstruct%edge_s
854 edge_n => gridstruct%edge_n
856 grid => gridstruct%grid
857 agrid => gridstruct%agrid
858 dxa => gridstruct%dxa
859 dya => gridstruct%dya
862 if (gridstruct%grid_type < 3)
then 872 if ( i==1 .and. j==1 )
goto 123
873 if ( i==2 .and. j==1 )
then 874 qin(0,-1) = qin(-1,2)
875 qin(0, 0) = qin(-1,1)
877 if ( i==1 .and. j==2 )
then 878 qin(-1,0) = qin(2,-1)
879 qin( 0,0) = qin(1,-1)
881 if ( i==2 .and. j==2 )
then 885 if ( i==npx .and. j==1 )
goto 123
886 if ( i==npx-1 .and. j==1 )
then 887 qin(npx,-1) = qin(npx+1,2)
888 qin(npx, 0) = qin(npx+1,1)
890 if ( i==npx-1 .and. j==2 )
then 891 qin(npx,0) = qin(npx-4,4)
893 if ( i==npx .and. j==2 )
then 894 qin(npx+1,0) = qin(npx-2,-1)
895 qin(npx, 0) = qin(npx-1,-1)
898 if ( i==npx .and. j==npy )
goto 123
899 if ( i==npx-1 .and. j==npy-1 )
then 900 qin(npx,npy) = qin(npx-4,npy-4)
902 if ( i==npx .and. j==npy-1 )
then 903 qin(npx+1,npy) = qin(npx-2,npy+1)
904 qin(npx, npy) = qin(npx-1,npy+1)
906 if ( i==npx-1 .and. j==npy )
then 907 qin(npx,npy+1) = qin(npx+1,npy-2)
908 qin(npx,npy ) = qin(npx+1,npy-1)
911 if ( i==1 .and. j==npy )
goto 123
912 if ( i==1 .and. j==npy-1 )
then 913 qin(-1,npy) = qin(2,npy+1)
914 qin( 0,npy) = qin(1,npy+1)
916 if ( i==2 .and. j==npy-1 )
then 917 qin(0,npy) = qin(4,npy-4)
919 if ( i==2 .and. j==npy )
then 920 qin(0,npy+1) = qin(-1,npy-2)
921 qin(0,npy ) = qin(-1,npy-1)
924 qout(i,j) = van2(1, i,j)*qin(i-2,j-2) + van2(2, i,j)*qin(i-1,j-2) + &
925 van2(3, i,j)*qin(i ,j-2) + van2(4, i,j)*qin(i+1,j-2) + &
926 van2(5, i,j)*qin(i-2,j-1) + van2(6, i,j)*qin(i-1,j-1) + &
927 van2(7, i,j)*qin(i ,j-1) + van2(8, i,j)*qin(i+1,j-1) + &
928 van2(9, i,j)*qin(i-2,j ) + van2(10,i,j)*qin(i-1,j ) + &
929 van2(11,i,j)*qin(i ,j ) + van2(12,i,j)*qin(i+1,j ) + &
930 van2(13,i,j)*qin(i-2,j+1) + van2(14,i,j)*qin(i-1,j+1) + &
931 van2(15,i,j)*qin(i ,j+1) + van2(16,i,j)*qin(i+1,j+1)
937 if ( gridstruct%sw_corner )
then 938 p0(1:2) = grid(1,1,1:2)
939 qout(1,1) = (
extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
940 extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
941 extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*
r3 944 if ( gridstruct%se_corner )
then 945 p0(1:2) = grid(npx,1,1:2)
946 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)) + &
947 extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
948 extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*
r3 950 if ( gridstruct%ne_corner )
then 951 p0(1:2) = grid(npx,npy,1:2)
952 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)) + &
953 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)) + &
954 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 956 if ( gridstruct%nw_corner )
then 957 p0(1:2) = grid(1,npy,1:2)
958 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)) + &
959 extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
960 extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*
r3 971 qx(i,j) =
b1*(qin(i-1,j)+qin(i,j)) +
b2*(qin(i-2,j)+qin(i+1,j))
977 qy(i,j) =
b1*(qin(i,j-1)+qin(i,j)) +
b2*(qin(i,j-2)+qin(i,j+1))
983 qout(i,j) = 0.5*(
a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
984 a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
990 if (
present(replace) )
then real, parameter b1
0.58333333
The type 'fv_grid_type' is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
subroutine, public a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
integer, parameter, public r_grid
The module 'a2b_edge' performs FV-consistent interpolation of pressure to corners.
The module 'fv_arrays' contains the 'fv_atmos_type' and associated datatypes.
real function, public great_circle_dist(q1, q2, radius)
The module 'fv_grid_utils' contains routines for setting up and computing grid-related quantities...
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
real function extrap_corner(p0, p1, p2, q1, q2)