Line data Source code
1 : !==================================================================================
2 : ! The subroutine reconstruction is called from both a horizontal
3 : ! threaded region and a nested region for horizontal and
4 : ! vertical threading. if horz_num_threads != horz_num_threads*vert_num_threads
5 : ! then the timers calls will generate segfault... So the simple solution is
6 : ! to deactivate them by default.
7 : !
8 : #define FVM_TIMERS .FALSE.
9 : !==================================================================================
10 : !MODULE FVM_RECONSTRUCTION_MOD--------------------------------------CE-for FVM!
11 : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011 !
12 : ! This module contains everything to do (ONLY) a CUBIC (3rd order) reconstruction !
13 : ! !
14 : ! IMPORTANT: the implementation is done for a ncfl > 1, which is not working !
15 : ! but it works for ncfl=1 !
16 : !
17 : ! This module has been recoded for multi-tracer efficiency (May, 2014)
18 : !
19 : !---------------------------------------------------------------------------!
20 : module fvm_reconstruction_mod
21 :
22 : use shr_kind_mod, only: r8=>shr_kind_r8
23 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
24 : use cam_abortutils, only: endrun
25 : use perf_mod, only: t_startf, t_stopf
26 :
27 :
28 : implicit none
29 : private
30 : ! integer, parameter, private:: nh = nhr+(nhe-1) ! = 2 (nhr=2; nhe=1)
31 : ! = 3 (nhr=2; nhe=2)
32 : public :: reconstruction, recons_val_cart, extend_panel_interpolate
33 : !reconstruction_gradient,
34 : contains
35 : ! ----------------------------------------------------------------------------------!
36 : !SUBROUTINE RECONSTRUCTION------------------------------------------------CE-for FVM!
37 : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011 !
38 : ! DESCRIPTION: controls the cubic (3rd order) reconstructions: !
39 : ! !
40 : ! CALLS: fillhalo_cubic, reconstruction_cubic !
41 : ! INPUT: fcube ... tracer values incl. the halo zone !
42 : ! fvm ... structure incl. tracer values aso ! !
43 : ! OUTPUT:recons ... has the reconstruction coefficients (5) for the 3rd order !
44 : ! reconstruction: dx, dy, dx^2, dy^2, dxdy !
45 : !-----------------------------------------------------------------------------------!
46 483116400 : subroutine reconstruction(fcube,nlev_in,k_in,recons,irecons,llimiter,ntrac_in,&
47 : nc,nhe,nhr,nhc,nht,ns,nh,&
48 : jx_min,jx_max,jy_min,jy_max,&
49 483116400 : cubeboundary,halo_interp_weight,ibase,&
50 483116400 : spherecentroid,&
51 483116400 : recons_metrics,recons_metrics_integral,&
52 483116400 : rot_matrix,centroid_stretch,&
53 483116400 : vertex_recons_weights,vtx_cart,&
54 : irecons_actual_in)
55 : implicit none
56 : !
57 : ! dimension(1-nhc:nc+nhc, 1-nhc:nc+nhc)
58 : !
59 : integer, intent(in) :: irecons
60 : integer, intent(in) :: nlev_in, k_in
61 : integer, intent(in) :: ntrac_in,nc,nhe,nhr,nhc,nht,ns,nh,cubeboundary
62 : real (kind=r8), dimension(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev_in,ntrac_in), intent(inout) :: fcube
63 : real (kind=r8), dimension(irecons,1-nhe:nc+nhe,1-nhe:nc+nhe,ntrac_in), intent(out) :: recons
64 : integer, intent(in) :: jx_min(3), jx_max(3), jy_min(3), jy_max(3)
65 : integer , intent(in):: ibase(1-nh:nc+nh,1:nhr,2)
66 : real (kind=r8), intent(in):: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
67 : real (kind=r8), intent(in):: spherecentroid(irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe)
68 : real (kind=r8), intent(in):: recons_metrics(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
69 : real (kind=r8), intent(in):: recons_metrics_integral(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
70 : integer , intent(in):: rot_matrix(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
71 : real (kind=r8), intent(in):: centroid_stretch(7,1-nhe:nc+nhe,1-nhe:nc+nhe)
72 : real (kind=r8), intent(in):: vertex_recons_weights(4,1:irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe)
73 : real (kind=r8), intent(in):: vtx_cart(4,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
74 :
75 : logical, intent(in) :: llimiter(ntrac_in)
76 : integer, optional, intent(in) :: irecons_actual_in
77 :
78 : integer :: irecons_actual
79 :
80 966232800 : real (kind=r8), dimension(1-nht:nc+nht,1-nht:nc+nht,3) :: f
81 :
82 : integer :: i,j,in,h,itr
83 : integer, dimension(2,3) :: jx,jy
84 :
85 483116400 : if (present(irecons_actual_in)) then
86 483116400 : irecons_actual = irecons_actual_in
87 : else
88 0 : irecons_actual = irecons
89 : end if
90 :
91 :
92 483116400 : jx(1,1)=jx_min(1); jx(2,1)=jx_max(1)-1
93 483116400 : jx(1,2)=jx_min(2); jx(2,2)=jx_max(2)-1
94 483116400 : jx(1,3)=jx_min(3); jx(2,3)=jx_max(3)-1
95 :
96 483116400 : jy(1,1)=jy_min(1); jy(2,1)=jy_max(1)-1
97 483116400 : jy(1,2)=jy_min(2); jy(2,2)=jy_max(2)-1
98 483116400 : jy(1,3)=jy_min(3); jy(2,3)=jy_max(3)-1
99 :
100 : !
101 : ! Initialize recons
102 : !
103 483116400 : call zero_non_existent_ghost_cell(recons,irecons,cubeboundary,nc,nhe,ntrac_in)
104 483116400 : if (irecons_actual>1) then
105 : if(FVM_TIMERS) call t_startf('FVM:reconstruction:part#1')
106 483116400 : if (nhe>0) then
107 5797396800 : do itr=1,ntrac_in
108 : call extend_panel_interpolate(nc,nhc,nhr,nht,ns,nh,&
109 5314280400 : fcube(:,:,k_in,itr),cubeboundary,halo_interp_weight,ibase,f(:,:,1),f(:,:,2:3))
110 : call get_gradients(f(:,:,:),jx,jy,irecons,recons(:,:,:,itr),&
111 5797396800 : rot_matrix,centroid_stretch,nc,nht,nhe,nhc,irecons_actual)
112 : end do
113 : else
114 0 : do itr=1,ntrac_in
115 : call extend_panel_interpolate(nc,nhc,nhr,nht,ns,nh,&
116 0 : fcube(:,:,k_in,itr),cubeboundary,halo_interp_weight,ibase,f(:,:,1))
117 : call get_gradients(f(:,:,:),jx,jy,irecons,recons(:,:,:,itr),&
118 0 : rot_matrix,centroid_stretch,nc,nht,nhe,nhc,irecons_actual)
119 : end do
120 : end if
121 : if(FVM_TIMERS) call t_stopf('FVM:reconstruction:part#1')
122 : if(FVM_TIMERS) call t_startf('FVM:reconstruction:part#2')
123 : !
124 : ! fill in non-existent (in physical space) corner values to simplify
125 : ! logic in limiter code (min/max operation)
126 : !
127 5797396800 : do itr=1,ntrac_in
128 5797396800 : if (llimiter(itr)) then
129 5314280400 : if (cubeboundary>4) then
130 5904756 : select case(cubeboundary)
131 : case (nwest)
132 17714268 : do h=1,nhe+1
133 11809512 : fcube(0,nc+h ,k_in,itr) = fcube(1-h,nc ,k_in,itr)
134 17714268 : fcube(1-h,nc+1,k_in,itr) = fcube(1 ,nc+h,k_in,itr)
135 : end do
136 : case (swest)
137 17714268 : do h=1,nhe+1
138 11809512 : fcube(1-h,0,k_in,itr) = fcube(1,1-h,k_in,itr)
139 17714268 : fcube(0,1-h,k_in,itr) = fcube(1-h,1,k_in,itr)
140 : end do
141 : case (seast)
142 17714268 : do h=1,nhe+1
143 11809512 : fcube(nc+h,0 ,k_in,itr) = fcube(nc,1-h,k_in,itr)
144 17714268 : fcube(nc+1,1-h,k_in,itr) = fcube(nc+h,1,k_in,itr)
145 : end do
146 : case (neast)
147 35428536 : do h=1,nhe+1
148 11809512 : fcube(nc+h,nc+1,k_in,itr) = fcube(nc,nc+h,k_in,itr)
149 17714268 : fcube(nc+1,nc+h,k_in,itr) = fcube(nc+h,nc,k_in,itr)
150 : end do
151 : end select
152 : end if
153 : call slope_limiter(nhe,nc,nhc,fcube(:,:,k_in,itr),jx,jy,irecons,recons(:,:,:,itr),&
154 : spherecentroid(:,1-nhe:nc+nhe,1-nhe:nc+nhe),&
155 5314280400 : recons_metrics,vertex_recons_weights,vtx_cart,irecons_actual)
156 : end if
157 : end do
158 : if(FVM_TIMERS) call t_stopf('FVM:reconstruction:part#2')
159 : end if
160 : if(FVM_TIMERS) call t_startf('FVM:reconstruction:part#3')
161 : select case (irecons_actual)
162 : case(1)
163 0 : do in=1,3
164 0 : do j=jy(1,in),jy(2,in)
165 0 : do i=jx(1,in),jx(2,in)
166 0 : recons(1,i,j,1:ntrac_in) = fcube(i,j,k_in,1:ntrac_in)
167 0 : recons(2:irecons,i,j,1:ntrac_in) = 0.0_r8
168 : end do
169 : end do
170 : end do
171 : case(3)
172 : ! do j=1-nhe,nc+nhe
173 : ! do i=1-nhe,nc+nhe
174 0 : do in=1,3
175 0 : do j=jy(1,in),jy(2,in)
176 0 : do i=jx(1,in),jx(2,in)
177 0 : recons(1,i,j,1:ntrac_in) = fcube(i,j,k_in,1:ntrac_in) &
178 0 : - recons(2,i,j,1:ntrac_in)*spherecentroid(1,i,j) &
179 0 : - recons(3,i,j,1:ntrac_in)*spherecentroid(2,i,j)
180 0 : recons(2,i,j,1:ntrac_in) = recons(2,i,j,1:ntrac_in)
181 0 : recons(3,i,j,1:ntrac_in) = recons(3,i,j,1:ntrac_in)
182 0 : recons(4:irecons,i,j,1:ntrac_in) = 0.0_r8
183 : end do
184 : end do
185 : end do
186 : case(6)
187 5797396800 : do itr=1,ntrac_in
188 : ! do j=1-nhe,nc+nhe
189 : ! do i=1-nhe,nc+nhe
190 21740238000 : do in=1,3
191 49576331376 : do j=jy(1,in),jy(2,in)
192 >17709*10^7 : do i=jx(1,in),jx(2,in)
193 >53133*10^7 : recons(1,i,j,itr) = fcube(i,j,k_in,itr) &
194 >13283*10^7 : - recons(2,i,j,itr)*spherecentroid(1,i,j) &
195 >13283*10^7 : - recons(3,i,j,itr)*spherecentroid(2,i,j) &
196 >13283*10^7 : + recons(4,i,j,itr)*recons_metrics_integral(1,i,j) &
197 >13283*10^7 : + recons(5,i,j,itr)*recons_metrics_integral(2,i,j) &
198 >66416*10^7 : + recons(6,i,j,itr)*recons_metrics_integral(3,i,j)
199 : recons(2,i,j,itr) = recons(2,i,j,itr) &
200 : - recons(4,i,j,itr)*2.0_r8*spherecentroid(1,i,j) &
201 >13283*10^7 : - recons(6,i,j,itr) *spherecentroid(2,i,j)
202 : recons(3,i,j,itr) = recons(3,i,j,itr) &
203 : - recons(5,i,j,itr)*2.0_r8*spherecentroid(2,i,j) &
204 >16115*10^7 : - recons(6,i,j,itr)*spherecentroid(1,i,j)
205 : !
206 : ! recons(i,j,4:6) already set in get_gradients
207 : !
208 : end do
209 : end do
210 : end do
211 : end do
212 : case default
213 483116400 : write(*,*) "irecons out of range in get_ceof", irecons
214 : end select
215 : if(FVM_TIMERS) call t_stopf('FVM:reconstruction:part#3')
216 :
217 : ! recons(a,b,3) * (centroid(a,b,1)**2 - centroid(a,b,3)) + &
218 : ! recons(a,b,4) * (centroid(a,b,2)**2 - centroid(a,b,4)) + &
219 : ! recons(a,b,5) * (centroid(a,b,1) * centroid(a,b,2) - centroid(a,b,5)) + &
220 :
221 :
222 : ! call debug_halo(fvm,fcubenew,fpanel)
223 : ! call debug_halo_recons(fvm,recons,recons_trunk)
224 : ! call print_which_case(fvm)
225 : !
226 : ! call debug_halo_neighbor (fvm,fotherface,fotherpanel)
227 : ! call debug_halo_neighbor_recons(fvm,recons,recons_trunk)
228 483116400 : end subroutine reconstruction
229 : !END SUBROUTINE RECONSTRUCTION--------------------------------------------CE-for FVM!
230 :
231 5314280400 : subroutine get_gradients(f,jx,jy,irecons,gradient,rot_matrix,centroid_stretch,nc,nht,nhe,nhc,irecons_actual)
232 : implicit none
233 : integer, intent(in) :: irecons,nc,nht,nhe,nhc,irecons_actual
234 : real (kind=r8), dimension(1-nht:nc+nht,1-nht:nc+nht,3), intent(in) :: f
235 : real (kind=r8), dimension(irecons,1-nhe:nc+nhe,1-nhe:nc+nhe), intent(inout):: gradient
236 : integer, dimension(2,3), intent(in) :: jx,jy
237 : integer , dimension(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc), intent(in) :: rot_matrix
238 : real (kind=r8), dimension(7,1-nhe:nc+nhe,1-nhe:nc+nhe), intent(in) :: centroid_stretch
239 :
240 : integer :: i,j,in
241 : real (kind=r8), dimension(2):: g
242 : real (kind=r8) :: sign
243 : character(len=128) :: errormsg
244 :
245 :
246 5314280400 : select case (irecons_actual)
247 : case(3)
248 0 : in=1
249 0 : do j=jy(1,in),jy(2,in)
250 0 : do i=jx(1,in),jx(2,in)
251 : !
252 : ! df/dx: 4-th-order finite difference: (-f(i+2)+8f(i+1)-8f(i-1)+f(i-2))/12dx
253 : !
254 0 : gradient(2,i,j) = -f(i+2,j ,in)+8.0_r8*f(i+1,j ,in)-8.0_r8*f(i-1,j ,in)+f(i-2,j ,in)
255 0 : gradient(3,i,j) = -f(i ,j+2,in)+8.0_r8*f(i ,j+1,in)-8.0_r8*f(i ,j-1,in)+f(i ,j-2,in)
256 : end do
257 : end do
258 0 : do in=2,3
259 0 : do j=jy(1,in),jy(2,in)
260 0 : do i=jx(1,in),jx(2,in)
261 0 : g(1) = -f(i+2,j ,in)+8.0_r8*f(i+1,j ,in)-8.0_r8*f(i-1,j ,in)+f(i-2,j ,in)
262 0 : g(2) = -f(i ,j+2,in)+8.0_r8*f(i ,j+1,in)-8.0_r8*f(i ,j-1,in)+f(i ,j-2,in)
263 0 : gradient(2:3,i,j) = MATMUL(rot_matrix(:,:,i,j),g(:))
264 : end do
265 : end do
266 : end do
267 0 : gradient(2,:,:) = centroid_stretch(1,:,:)*gradient(2,:,:)
268 0 : gradient(3,:,:) = centroid_stretch(2,:,:)*gradient(3,:,:)
269 : case (6)
270 5314280400 : in=1
271 31531397040 : do j=jy(1,in),jy(2,in)
272 >16086*10^7 : do i=jx(1,in),jx(2,in)
273 : !
274 : ! df/dx: 4-th-order finite difference: (-f(i+2)+8f(i+1)-8f(i-1)+f(i-2))/12dx
275 : !
276 >12933*10^7 : gradient(2,i,j) = -f(i+2,j ,in)+ 8.0_r8*f(i+1,j ,in) - 8.0_r8*f(i-1,j ,in)+f(i-2,j ,in)
277 >12933*10^7 : gradient(3,i,j) = -f(i ,j+2,in)+ 8.0_r8*f(i ,j+1,in) - 8.0_r8*f(i ,j-1,in)+f(i ,j-2,in)
278 : !
279 : ! d2f/dx2:
280 : !
281 >12933*10^7 : gradient(4,i,j) = -f(i+2,j ,in)+16.0_r8*f(i+1,j ,in)-30.0_r8*f(i,j,in)+16.0_r8*f(i-1,j ,in)-f(i-2,j ,in)
282 >12933*10^7 : gradient(5,i,j) = -f(i ,j+2,in)+16.0_r8*f(i ,j+1,in)-30.0_r8*f(i,j,in)+16.0_r8*f(i ,j-1,in)-f(i ,j-2,in)
283 :
284 >12933*10^7 : gradient(6,i,j) = f(i+1,j+1,in)- f(i+1,j-1,in) - f(i-1,j+1,in)+f(i-1,j-1,in)
285 : !
286 : ! "stretching factors
287 : !
288 >12933*10^7 : gradient(2,i,j) = centroid_stretch(1,i,j)*gradient(2,i,j)
289 >12933*10^7 : gradient(3,i,j) = centroid_stretch(2,i,j)*gradient(3,i,j)
290 :
291 >12933*10^7 : gradient(4,i,j) = centroid_stretch(3,i,j)*gradient(4,i,j)+centroid_stretch(6,i,j)*gradient(2,i,j)
292 >12933*10^7 : gradient(5,i,j) = centroid_stretch(4,i,j)*gradient(5,i,j)+centroid_stretch(7,i,j)*gradient(3,i,j)
293 :
294 >15555*10^7 : gradient(6,i,j) = centroid_stretch(5,i,j)*gradient(6,i,j)
295 : end do
296 : end do
297 15942841200 : do in=2,3
298 74399925600 : if (SUM(rot_matrix(:,:,jx(1,in),jy(1,in)))==0) then
299 : sign=-1
300 : else
301 10282148448 : sign=1
302 : end if
303 18044934336 : do j=jy(1,in),jy(2,in)
304 16226269488 : do i=jx(1,in),jx(2,in)
305 3495615552 : g(1) = -f(i+2,j ,in)+8.0_r8*f(i+1,j ,in)-8.0_r8*f(i-1,j ,in)+f(i-2,j ,in)
306 3495615552 : g(2) = -f(i ,j+2,in)+8.0_r8*f(i ,j+1,in)-8.0_r8*f(i ,j-1,in)+f(i ,j-2,in)
307 24469308864 : gradient(2:3,i,j) = MATMUL(rot_matrix(:,:,i,j),g(:))
308 :
309 3495615552 : g(1) = -f(i+2,j ,in)+16.0_r8*f(i+1,j ,in)-30.0_r8*f(i,j,in)+16.0_r8*f(i-1,j ,in)-f(i-2,j ,in)
310 3495615552 : g(2) = -f(i ,j+2,in)+16.0_r8*f(i ,j+1,in)-30.0_r8*f(i,j,in)+16.0_r8*f(i ,j-1,in)-f(i ,j-2,in)
311 24469308864 : gradient(4:5,i,j) = MATMUL(ABS(rot_matrix(:,:,i,j)),g(:))
312 :
313 3495615552 : gradient(6,i,j) = sign*(f(i+1,j+1,in)- f(i+1,j-1,in) - f(i-1,j+1,in)+f(i-1,j-1,in))
314 : !
315 : ! "stretching factors
316 : !
317 3495615552 : gradient(2,i,j) = centroid_stretch(1,i,j)*gradient(2,i,j)
318 3495615552 : gradient(3,i,j) = centroid_stretch(2,i,j)*gradient(3,i,j)
319 :
320 3495615552 : gradient(4,i,j) = centroid_stretch(3,i,j)*gradient(4,i,j)+centroid_stretch(6,i,j)*gradient(2,i,j)
321 3495615552 : gradient(5,i,j) = centroid_stretch(4,i,j)*gradient(5,i,j)+centroid_stretch(7,i,j)*gradient(3,i,j)
322 :
323 5597708688 : gradient(6,i,j) = centroid_stretch(5,i,j)*gradient(6,i,j)
324 : end do
325 : end do
326 : end do
327 : case default
328 0 : write(errormsg, *) irecons
329 5314280400 : call endrun('ERROR: irecons out of range in slope_limiter'//errormsg)
330 : end select
331 5314280400 : end subroutine get_gradients
332 :
333 :
334 5314280400 : subroutine slope_limiter(nhe,nc,nhc,fcube,jx,jy,irecons,recons,spherecentroid,recons_metrics,&
335 5314280400 : vertex_recons_weights,vtx_cart,irecons_actual)
336 : implicit none
337 : integer , intent(in) :: irecons,nhe,nc,nhc,irecons_actual
338 : real (kind=r8), dimension(1-nhc:, 1-nhc:) , intent(in) :: fcube
339 : real (kind=r8), dimension(irecons,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(inout):: recons
340 : integer, dimension(2,3) , intent(in) :: jx,jy
341 : real (kind=r8), dimension(irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(in) :: spherecentroid
342 : real (kind=r8), dimension(3,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(in) :: recons_metrics
343 : real (kind=r8), dimension(4,1:irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe), intent(in) :: vertex_recons_weights
344 : real (kind=r8), dimension(4,2,1-nhc:nc+nhc,1-nhc:nc+nhc) , intent(in) :: vtx_cart
345 :
346 : real (kind=r8):: minval_patch,maxval_patch
347 : real (kind=r8):: phi, min_val, max_val,disc
348 :
349 : real (kind=r8):: min_phi
350 : real (kind=r8):: extrema(2), xminmax(2),yminmax(2),extrema_value(13)
351 :
352 : real(kind=r8) :: invtmp ! temporary to pre-compute inverses
353 : integer :: itmp1,itmp2,i,j,in,vertex,n
354 :
355 : ! real (kind=r8), dimension(-1:5) :: diff_value
356 : real (kind=r8), dimension(-1:1) :: minval_array, maxval_array
357 : real (kind=r8), parameter :: threshold = 1.0E-40_r8
358 : character(len=128) :: errormsg
359 5314280400 : select case (irecons_actual)
360 : !
361 : ! PLM limiter
362 : !
363 : case(3)
364 0 : do in=1,3
365 0 : do j=jy(1,in),jy(2,in)
366 0 : do i=jx(1,in),jx(2,in)
367 : !rck combined min/max and unrolled inner loop
368 : !minval_patch = MINVAL(fcube(i-1:i+1,j-1:j+1))
369 : !maxval_patch = MAXVAL(fcube(i-1:i+1,j-1:j+1))
370 : !DIR$ SIMD
371 0 : do itmp2=-1,+1
372 0 : itmp1 = j+itmp2
373 0 : minval_array(itmp2) = min(fcube(i-1,itmp1),fcube(i,itmp1),fcube(i+1,itmp1))
374 0 : maxval_array(itmp2) = max(fcube(i-1,itmp1),fcube(i,itmp1),fcube(i+1,itmp1))
375 : enddo
376 0 : minval_patch = min(minval_array(-1),minval_array(0),minval_array(1))
377 0 : maxval_patch = max(maxval_array(-1),maxval_array(0),maxval_array(1))
378 :
379 0 : min_phi=1.0_r8
380 :
381 : !
382 : ! coordinate bounds (could be pre-computed!)
383 : !
384 0 : xminmax(1) = min(vtx_cart(1,1,i,j),vtx_cart(2,1,i,j),vtx_cart(3,1,i,j),vtx_cart(4,1,i,j))
385 0 : xminmax(2) = max(vtx_cart(1,1,i,j),vtx_cart(2,1,i,j),vtx_cart(3,1,i,j),vtx_cart(4,1,i,j))
386 0 : yminmax(1) = min(vtx_cart(1,2,i,j),vtx_cart(2,2,i,j),vtx_cart(3,2,i,j),vtx_cart(4,2,i,j))
387 0 : yminmax(2) = max(vtx_cart(1,2,i,j),vtx_cart(2,2,i,j),vtx_cart(3,2,i,j),vtx_cart(4,2,i,j))
388 :
389 : !rck restructured loop
390 : !DIR$ SIMD
391 0 : do vertex=1,4
392 0 : call recons_val_cart_plm(fcube(i,j), vtx_cart(vertex,1,i,j), vtx_cart(vertex,2,i,j), spherecentroid(:,i,j), &
393 0 : recons(1:3,i,j), extrema_value(vertex))
394 : end do
395 0 : max_val = MAXVAL(extrema_value(1:4))
396 0 : min_val = MINVAL(extrema_value(1:4))
397 :
398 0 : if (max_val>maxval_patch.and.abs(max_val-fcube(i,j))>threshold) then
399 0 : phi = (maxval_patch-fcube(i,j))/(max_val-fcube(i,j))
400 0 : if (phi<min_phi) min_phi=phi
401 : end if
402 0 : if (min_val<minval_patch.and.abs(min_val-fcube(i,j))>threshold) then
403 0 : phi = (minval_patch-fcube(i,j))/(min_val-fcube(i,j))
404 0 : if (phi<min_phi) min_phi=phi
405 : end if
406 : ! Apply monotone limiter to all reconstruction coefficients
407 0 : recons(2:3,i,j)=min_phi*recons(2:3,i,j)
408 : end do
409 : end do
410 : end do
411 : !
412 : ! PPM limiter
413 : !
414 : case(6)
415 : !
416 : ! default branch
417 : !
418 21257121600 : do in=1,3
419 49576331376 : do j=jy(1,in),jy(2,in)
420 >17709*10^7 : do i=jx(1,in),jx(2,in)
421 : !rck combined min/max and unrolled inner loop
422 : !minval_patch = MINVAL(fcube(i-1:i+1,j-1:j+1))
423 : !maxval_patch = MAXVAL(fcube(i-1:i+1,j-1:j+1))
424 : !DIR$ SIMD
425 >53133*10^7 : do itmp2=-1,+1
426 >39850*10^7 : itmp1 = j+itmp2
427 >39850*10^7 : minval_array(itmp2) = min(fcube(i-1,itmp1),fcube(i,itmp1),fcube(i+1,itmp1))
428 >53133*10^7 : maxval_array(itmp2) = max(fcube(i-1,itmp1),fcube(i,itmp1),fcube(i+1,itmp1))
429 : enddo
430 >13283*10^7 : minval_patch = min(minval_array(-1),minval_array(0),minval_array(1))
431 >13283*10^7 : maxval_patch = max(maxval_array(-1),maxval_array(0),maxval_array(1))
432 :
433 >13283*10^7 : min_phi=1.0_r8
434 : !rck restructured loop
435 :
436 >66416*10^7 : extrema_value(1:4) = fcube(i,j)
437 : !DIR$ SIMD
438 >66416*10^7 : do vertex=1,4
439 >33208*10^8 : do itmp1=1,irecons-1
440 >26566*10^8 : extrema_value(vertex) = extrema_value(vertex) + &
441 >58446*10^8 : recons(itmp1+1,i,j)*vertex_recons_weights(vertex,itmp1,i,j)
442 : enddo
443 : enddo
444 >13283*10^8 : extrema_value(5:13) = extrema_value(1)
445 : !
446 : ! coordinate bounds (could be pre-computed!)
447 : !
448 >13283*10^7 : xminmax(1) = min(vtx_cart(1,1,i,j),vtx_cart(2,1,i,j),vtx_cart(3,1,i,j),vtx_cart(4,1,i,j))
449 >13283*10^7 : xminmax(2) = max(vtx_cart(1,1,i,j),vtx_cart(2,1,i,j),vtx_cart(3,1,i,j),vtx_cart(4,1,i,j))
450 >13283*10^7 : yminmax(1) = min(vtx_cart(1,2,i,j),vtx_cart(2,2,i,j),vtx_cart(3,2,i,j),vtx_cart(4,2,i,j))
451 >13283*10^7 : yminmax(2) = max(vtx_cart(1,2,i,j),vtx_cart(2,2,i,j),vtx_cart(3,2,i,j),vtx_cart(4,2,i,j))
452 :
453 : ! Check if the quadratic is minimized within the element
454 : ! Extrema in the interior of the element (there might be just one candidate)
455 : ! DO NOT NEED ABS here, if disc<0 we have a saddle point (no maximum or minimum)
456 >13283*10^7 : disc = 4.0_r8 * recons(4,i,j) * recons(5,i,j) - recons(6,i,j)**2
457 >13283*10^7 : if (abs(disc) > threshold) then
458 78103680914 : extrema(1) = recons(6,i,j) * recons(3,i,j) - 2.0_r8 * recons(5,i,j) * recons(2,i,j)
459 78103680914 : extrema(2) = recons(6,i,j) * recons(2,i,j) - 2.0_r8 * recons(4,i,j) * recons(3,i,j)
460 :
461 78103680914 : disc=1.0_r8/disc
462 78103680914 : extrema(1) = extrema(1) * disc + spherecentroid(1,i,j)
463 78103680914 : extrema(2) = extrema(2) * disc + spherecentroid(2,i,j)
464 : if ( (extrema(1) - xminmax(1) > -threshold) .and. & !xmin
465 : (extrema(1) - xminmax(2) < threshold) .and. & !xmax
466 78103680914 : (extrema(2) - yminmax(1) > -threshold) .and. & !ymin
467 : (extrema(2) - yminmax(2) < threshold)) then !ymax
468 : call recons_val_cart(fcube(i,j), extrema(1), extrema(2), spherecentroid(:,i,j), &
469 11583977935 : recons_metrics(:,i,j), recons(:,i,j), extrema_value(5))
470 : endif
471 : endif
472 : !
473 : ! Check all potential minimizer points along element boundaries
474 : !
475 >13283*10^7 : if (abs(recons(6,i,j)) > threshold) then
476 >10316*10^7 : invtmp = 1.0_r8 / (recons(6,i,j) + spherecentroid(2,i,j))
477 >30950*10^7 : do n=1,2
478 : ! Left edge, intercept with du/dx = 0
479 >20633*10^7 : extrema(2) = invtmp * (-recons(2,i,j) - 2.0_r8 * recons(4,i,j) * (xminmax(n) - spherecentroid(1,i,j)))
480 >30950*10^7 : if ((extrema(2) > yminmax(1)-threshold) .and. (extrema(2) < yminmax(2)+threshold)) then
481 : call recons_val_cart(fcube(i,j), xminmax(n), extrema(2), spherecentroid(:,i,j), &
482 2013757488 : recons_metrics(:,i,j), recons(:,i,j), extrema_value(5+n))
483 : endif
484 : enddo
485 : ! Top/bottom edge, intercept with du/dy = 0
486 >10316*10^7 : invtmp = 1.0_r8 / recons(6,i,j) + spherecentroid(1,i,j)
487 >30950*10^7 : do n = 1,2
488 >20633*10^7 : extrema(1) = invtmp * (-recons(3,i,j) - 2.0_r8 * recons(5,i,j) * (yminmax(n) - spherecentroid(2,i,j)))
489 >30950*10^7 : if ((extrema(1) > xminmax(1)-threshold) .and. (extrema(1) < xminmax(2)+threshold)) then
490 : call recons_val_cart(fcube(i,j), extrema(1), yminmax(n),spherecentroid(:,i,j), &
491 2023413045 : recons_metrics(:,i,j), recons(:,i,j), extrema_value(7+n))
492 : endif
493 : enddo
494 : endif
495 :
496 : ! Top/bottom edge, y=const., du/dx=0
497 >13283*10^7 : if (abs(recons(4,i,j)) > threshold) then
498 >10299*10^7 : invtmp = 1.0_r8 / (2.0_r8 * recons(4,i,j))! + spherecentroid(1,i,j)
499 >30897*10^7 : do n = 1,2
500 >20598*10^7 : extrema(1) = spherecentroid(1,i,j)+&
501 >20598*10^7 : invtmp * (-recons(2,i,j) - recons(6,i,j) * (yminmax(n) - spherecentroid(2,i,j)))
502 :
503 >30897*10^7 : if ((extrema(1) > xminmax(1)-threshold) .and. (extrema(1) < xminmax(2)+threshold)) then
504 : call recons_val_cart(fcube(i,j), extrema(1), yminmax(n), spherecentroid(:,i,j),&
505 44139495920 : recons_metrics(:,i,j),recons(:,i,j), extrema_value(9+n))
506 : endif
507 : enddo
508 : endif
509 : ! Left/right edge, x=const., du/dy=0
510 >13283*10^7 : if (abs(recons(5,i,j)) > threshold) then
511 >10323*10^7 : invtmp = 1.0_r8 / (2.0_r8 * recons(5,i,j))
512 >30971*10^7 : do n = 1,2
513 >20647*10^7 : extrema(2) = spherecentroid(2,i,j)+&
514 >20647*10^7 : invtmp * (-recons(3,i,j) - recons(6,i,j) * (xminmax(n) - spherecentroid(1,i,j)))
515 :
516 >30971*10^7 : if ((extrema(2)>yminmax(1)-threshold) .and. (extrema(2) < yminmax(2)+threshold)) then
517 : call recons_val_cart(fcube(i,j), xminmax(n), extrema(2), spherecentroid(:,i,j), &
518 51376500713 : recons_metrics(:,i,j), recons(:,i,j), extrema_value(11+n))
519 : endif
520 : enddo
521 : endif
522 : !rck - combined min/max calculation and unrolled
523 : ! max_val = MAXVAL(extrema_value)
524 : ! min_val = MINVAL(extrema_value)
525 >13283*10^7 : max_val = extrema_value(13)
526 >13283*10^7 : min_val = extrema_value(13)
527 >13283*10^7 : do itmp1 = 1,12,4
528 >39850*10^7 : max_val = max(max_val, extrema_value(itmp1),extrema_value(itmp1+1),extrema_value(itmp1+2),extrema_value(itmp1+3))
529 >39850*10^7 : min_val = min(min_val, extrema_value(itmp1),extrema_value(itmp1+1),extrema_value(itmp1+2),extrema_value(itmp1+3))
530 : enddo
531 : !rck
532 :
533 >13283*10^7 : if (max_val>maxval_patch.and.abs(max_val-fcube(i,j))>threshold) then
534 18311089416 : phi = (maxval_patch-fcube(i,j))/(max_val-fcube(i,j))
535 18311089416 : if (phi<min_phi) min_phi=phi
536 : end if
537 >13283*10^7 : if (min_val<minval_patch.and.abs(min_val-fcube(i,j))>threshold) then
538 35306905308 : phi = (minval_patch-fcube(i,j))/(min_val-fcube(i,j))
539 35306905308 : if (phi<min_phi) min_phi=phi
540 : end if
541 >82531*10^7 : recons(2:6,i,j)=min_phi*recons(2:6,i,j)
542 : end do
543 : end do
544 : end do
545 : case default
546 0 : write(errormsg, *) irecons
547 5314280400 : call endrun('ERROR: irecons out of range in slope_limiter'//errormsg)
548 : end select
549 :
550 5314280400 : end subroutine slope_limiter
551 :
552 :
553 :
554 : ! ----------------------------------------------------------------------------------!
555 : !SUBROUTINE RECONS_VAL_CART-----------------------------------------------CE-for FVM!
556 : ! AUTHOR: CHRISTOPH ERATH, 30.November 2011 !
557 : ! DESCRIPTION: returns the value from the reconstruction (3rd order Taylor polynom) !
558 : ! at the point (cartx,carty) -> in cube CARTESIAN coordinates !
559 : ! !
560 : ! INPUT: fcube ... tracer values incl. the halo zone !
561 : ! cartx ... x cartesian coordinate of the evaluation point !
562 : ! carty ... y cartesian coordinate of the evaluation point !
563 : ! centroid.. x,y,x^2,y^2,xy !
564 : ! recons ... array of reconstructed coefficients !
565 : ! OUTPUT: value ... evaluation at a given point !
566 : !-----------------------------------------------------------------------------------!
567 >11113*10^7 : subroutine recons_val_cart(fcube, cartx, carty, centroid, pre_computed_metrics, recons, value)
568 : implicit none
569 : real(kind=r8), intent(in) :: fcube
570 : real(kind=r8), intent(in) :: cartx, carty
571 : real(kind=r8), dimension(1:5), intent(in) :: centroid
572 : real(kind=r8), dimension(3), intent(in) :: pre_computed_metrics
573 : real(kind=r8), dimension(1:6), intent(in) :: recons
574 : real(kind=r8), intent(out) :: value
575 : real(kind=r8) :: dx, dy
576 >11113*10^7 : dx = cartx - centroid(1)
577 >11113*10^7 : dy = carty - centroid(2)
578 : ! Evaluate constant order terms
579 : value = fcube + &
580 : ! Evaluate linear order terms
581 : recons(2) * dx + &
582 : recons(3) * dy + &
583 : ! Evaluate second order terms
584 : recons(4) * (pre_computed_metrics(1) + dx*dx) + &
585 : recons(5) * (pre_computed_metrics(2) + dy*dy) + &
586 >11113*10^7 : recons(6) * (pre_computed_metrics(3) + dx*dy)
587 >11113*10^7 : END subroutine recons_val_cart
588 :
589 0 : subroutine recons_val_cart_plm(fcube, cartx, carty, centroid, recons, value)
590 : implicit none
591 : real(kind=r8), intent(in) :: fcube
592 : real(kind=r8), intent(in) :: cartx, carty
593 : real(kind=r8), dimension(1:5), intent(in) :: centroid
594 : real(kind=r8), dimension(1:3), intent(in) :: recons
595 : real(kind=r8), intent(out) :: value
596 : real(kind=r8) :: dx, dy
597 0 : dx = cartx - centroid(1)
598 0 : dy = carty - centroid(2)
599 : ! Evaluate constant order terms
600 : value = fcube + &
601 : ! Evaluate linear order terms
602 : recons(2) * dx + &
603 0 : recons(3) * dy
604 0 : END subroutine recons_val_cart_plm
605 :
606 :
607 : ! ----------------------------------------------------------------------------------!
608 : !SUBROUTINE SLOPELIMITER_VAL----------------------------------------------CE-for FVM!
609 : ! AUTHOR: CHRISTOPH ERATH, 30.November 2011 !
610 : ! DESCRIPTION: returns the value from the reconstruction (3rd order Taylor polynom) !
611 : ! at the point (cartx,carty) -> in cube CARTESIAN coordinates !
612 : ! !
613 : ! INPUT: value ... point value (calculated here by recons_val_cart) !
614 : ! cell_value ... tracer value (in the cell center) of the cell !
615 : ! local_min ... minmal value in the patch !
616 : ! local_max ... maximal value in the patch !
617 : ! INPUT/OUTPUT: min_phi ... slope limiter, inout because we go through any possible !
618 : ! extrema on the cell !
619 : !-----------------------------------------------------------------------------------!
620 : subroutine slopelimiter_val(value, cell_value, local_min, local_max, min_phi)
621 : implicit none
622 : real (kind=r8), intent(in) :: value, cell_value
623 : real (kind=r8), intent(in) :: local_min, local_max
624 : real (kind=r8), intent(inout) :: min_phi
625 : real (kind=r8) :: phi
626 :
627 : ! Check against the minimum bound on the reconstruction
628 : if (value - cell_value > 1.0e-12_r8 * value) then
629 : phi = (local_max - cell_value) / (value - cell_value)
630 : if (phi < min_phi) then
631 : min_phi = phi
632 : endif
633 : ! Check against the maximum bound on the reconstruction
634 : elseif (value - cell_value < -1.0e-12_r8 * value) then
635 : phi = (local_min - cell_value) / (value - cell_value)
636 : if(phi < min_phi) then
637 : min_phi = phi
638 : endif
639 : endif
640 : end subroutine slopelimiter_val
641 : !END SUBROUTINE SLOPELIMITER_VAL------------------------------------------CE-for FVM!
642 :
643 17147411424 : function matmul_w(w,f,ns)
644 : implicit none
645 : real (kind=r8) :: matmul_w
646 : real (kind=r8),dimension(:), intent(in) :: w,f !dimension(ns)
647 : integer, intent(in) :: ns
648 : integer :: k
649 17147411424 : matmul_w = 0.0_r8
650 68589645696 : do k=1,ns
651 68589645696 : matmul_w = matmul_w+w(k)*f(k)
652 : end do
653 17147411424 : end function matmul_w
654 :
655 : ! special hard-coded version of the function where ns=3
656 : ! for performance optimization
657 : ! function matmul_w(w, f)
658 : ! IMPLICIT NONE
659 : ! REAL(KIND=r8), dimension(3), intent(in) :: w
660 : ! REAL(KIND=r8), dimension(3), intent(in) :: f
661 : ! REAL(KIND=r8) :: matmul_w
662 : ! matmul_w = w(1)*f(1) + w(2)*f(2) + w(3)*f(3)
663 : ! end function matmul_w
664 :
665 5314280400 : subroutine extend_panel_interpolate(nc,nhc,nhr,nht,ns,nh,fcube,cubeboundary,halo_interp_weight,ibase,&
666 5314280400 : fpanel,fotherpanel)
667 : implicit none
668 : integer, intent(in) :: cubeboundary,nc,nhr,nht,nh,nhc,ns
669 : real (kind=r8), &
670 : dimension(1-nhc:nc+nhc, 1-nhc:nc+nhc), intent(in) :: fcube
671 :
672 : real (kind=r8), intent(in) :: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
673 : integer , intent(in) :: ibase(1-nh:nc+nh,1:nhr,2)
674 :
675 : real (kind=r8) , dimension(1-nht:nc+nht, 1-nht:nc+nht ), intent(out) :: fpanel
676 : real (kind=r8), dimension(1-nht:nc+nht,1-nht:nc+nht,2), intent(out), optional :: fotherpanel
677 :
678 : integer :: i, halo,ibaseref
679 10628560800 : real (kind=r8), dimension(1:ns,1-nh:nc+nh,1:nhr) :: w
680 : !
681 : ! fpanel = 1.0E19 !dbg
682 : !
683 : !
684 : ! Stencil for reconstruction is:
685 : !
686 : ! ---------------------
687 : ! | | | i | | |
688 : ! ---------------------
689 : ! | | i | i | i | |
690 : ! ---------------------
691 : ! | i | i | R | i | i |
692 : ! ---------------------
693 : ! | | i | i | i | |
694 : ! ---------------------
695 : ! | | | i | | |
696 : ! ---------------------
697 : !
698 : ! where
699 : !
700 : ! "R" is cell for which we whish to do the reconstruction
701 : ! "i" is the stencil
702 : !
703 : !
704 : ! If one or more point in the stencil is on another panel(s) then we need to interpolate
705 : ! to a stencil that is an extension of the panel on which R is located
706 : ! (this is done using one dimensional cubic Lagrange interpolation along panel side)
707 : !
708 : ! Example: say that southern most "s" on Figure above is on another panels projection then the stencil becomes
709 : !
710 : !
711 : ! ---------------------------------
712 : ! | | | | | | i | | |
713 : ! ----------------|----------------
714 : ! | | | | | i | i | i | |
715 : ! ----------------|----------------
716 : ! | | | | i | i | R | i | i |
717 : ! ----------------|----------------
718 : ! | | | | | i | i | i | |
719 : ! ---------------------------------
720 : ! / / / / / S /S&i/ S / S /
721 : ! /---/---/---/---/---/---/---/---/
722 : ! / / / / / / / / /
723 : !/---/---/---/---/---/---/---/---/
724 : !
725 : !
726 : ! where "S" are the cell average values used for the cubic interpolation (located on the South panel)
727 : !
728 : !
729 5314280400 : if (cubeboundary==0) then
730 >42126*10^7 : fpanel(1-nht:nc+nht,1-nht:nc+nht)=fcube(1-nht:nc+nht,1-nht:nc+nht)
731 : else if (cubeboundary==west) then
732 : ! !
733 : ! ! Case shown below: nhr=2, nhe=1, nht=nhr+nhe
734 : ! ! (nhr = reconstruction width along x and y)
735 : ! ! (nhe = max. Courant number)
736 : ! !
737 : ! !
738 : ! Figure below shows the element in question ! In terms of data structure:
739 : ! (center element) and the surrounding elements !
740 : ! on the element in question's projection ! * "H" is on same panel average value
741 : ! ! * "w" is west panel average values that need
742 : ! Notation: "0" marks the element boundaries ! to be interpolated to main element
743 : ! ! projection
744 : ! Elements to the west are on a different projection ! * "i" is extra halo required by the cubic
745 : ! ! interpolation
746 : ! 0 !
747 : ! |0000 !
748 : ! | |00000 !
749 : ! |\--| |000000000000000000000000000000000000 ! -x---x---x---x---x---x---x---x---x---x---x---x
750 : ! | |\--| 0 | | | 0 | | | 0 ! | | | i | | | | | | | | |
751 : ! |\--| |\--0---------------0---------------0 ! -------------x---------------x---------------x
752 : ! | |\--| 0 | | | 0 | | | 0 ! | | i | i | H | H | H | H | H | | | |
753 : ! |\--| |\--0---------------0---------------0 ! -------------x---------------x---------------x
754 : ! 0 |\--| 0 | | | 0 | | | 0 ! | | i | w | H | H | H | H | H | H | | |
755 : ! |0000 |\--0---------------0---------------0 ! -------------x---------------x---------------x
756 : ! | |0000 0 | | | 0 | | | 0 ! | | w | w | r | r | r | r | r | H | H | |
757 : ! |\--| |000000000000000000000000000000000000 ! -x---x---x---00000000000000000---x---x---x---x
758 : ! | |\--| 0 | | | 0 | | | 0 ! | | w | w 0 r | r | r | r 0 r | H | H | |
759 : ! |\--| \---0---------------0---------------0 ! -------------0---------------0---------------x
760 : ! | |\--| 0 | | | 0 | | | 0 ! | | w | w 0 r | r | r | r 0 r | H | H | |
761 : ! |\--| \---0---------------0---------------0 ! -------------0---------------0---------------x
762 : ! 0 |\--| 0 | | | 0 | | | 0 ! | | w | w 0 r | r | r | r 0 r | H | H | |
763 : ! |0000 |\--0---------------0---------------0 ! -------------0---------------0---------------x
764 : ! | |0000 0 | | | 0 | | | 0 ! | | w | w 0 r | r | r | r 0 r | H | H | |
765 : ! |\--| |000000000000000000000000000000000000 ! -x---x---x---00000000000000000---x---x---x---x
766 : ! | |\--| 0 | | | 0 | | | 0 ! | | w | w | r | r | r | r | r | H | H | |
767 : ! |\--| |\--0---------------0---------------0 ! -------------x---------------x---------------x
768 : ! | |\--| 0 | | | 0 | | | 0 ! | | i | w | H | H | H | H | H | H | | |
769 : ! |\--| |\--0---------------0---------------0 ! -------------x---------------x---------------x
770 : ! 0 |\--| 0 | | | 0 | | | 0 ! | | i | i | H | H | H | H | H | | | |
771 : ! 0000 |\--0---------------0---------------0 ! -------------x---------------x---------------x
772 : ! 0000 0 | | | 0 | | | 0 ! | | | i | | | | | | | | |
773 : ! 000000000000000000000000000000000000 ! -x---x---x---x---x---x---x---x---x---x---x---x
774 : !
775 : !
776 : ! -2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
777 : !
778 : !
779 : !
780 : ! fill in values (incl. halo) that are on the "main" panels projection
781 : !
782 10581322752 : fpanel(1:nc+nht,1-nht:nc+nht)=fcube(1:nc+nht,1-nht:nc+nht)
783 : !
784 : ! fill in values that are on the west panels projection
785 : !
786 9754656912 : w = halo_interp_weight(:,:,:,1)
787 495999504 : do halo=1,nhr
788 2479997520 : do i=halo-nh,nc+nh-(halo-1)
789 1983998016 : ibaseref=ibase(i,halo,1)
790 : ! ibaseref = ibase(i,halo,1)
791 2314664352 : fpanel(1-halo ,i) = matmul_w(w(:,i,halo),fcube(1-halo ,ibaseref:ibaseref+ns-1),ns)
792 : end do
793 : end do
794 :
795 165333168 : if (present(fotherpanel)) then
796 : !
797 : ! fill in values that are on the west panels projection
798 : !
799 6117327216 : fotherpanel (1-nht:0,1-nht:nc+nht,1)=fcube(1-nht:0,1-nht:nc+nht)
800 : !
801 495999504 : do halo=1,nhr
802 2479997520 : do i=halo-nh,nc+nh-(halo-1)
803 1983998016 : ibaseref=ibase(i,halo,1)
804 : !
805 : ! Exploit symmetry in interpolation weights
806 : !
807 2314664352 : fotherpanel(halo,i,1) = matmul_w(w(:,i,halo),fcube(halo ,ibaseref:ibaseref+ns-1),ns)
808 : end do
809 : end do
810 : end if
811 : else if (cubeboundary==east) then
812 : !
813 : ! north part is on different panel
814 : !
815 : ! stencil
816 : !
817 : ! CN<1 case !
818 : ! !
819 : !
820 : !
821 : ! 0 !
822 : ! 0000| !
823 : ! 0000| | !
824 : ! 000000000000000000000000000000000000| |--/| ! x---x---x---x---x---x---x---x---x---x---x---x-
825 : ! 0 | | | 0 | | | 0 |--/| | ! | | | | | | | | | i | | |
826 : ! 0---------------0---------------0--/ |--/| ! x---------------x---------------x---x---x---x-
827 : ! 0 | | | 0 | | | 0 |--/| | ! | | | | H | H | H | H | H | i | i | |
828 : ! 0---------------0---------------0--/| |--/| ! x---------------x---------------x---x---x---x-
829 : ! 0 | | | 0 | | | 0 |--/| 0 ! | | | H | H | H | H | H | H | e | i | |
830 : ! 0---------------0---------------0--/| 0000| ! x---------------x---------------x---x---x---x-
831 : ! 0 | | | 0 | | | 0 0000| | ! | | H | H | r | r | r | r | r | e | e | |
832 : ! 000000000000000000000000000000000000| |--/| ! x---x---x---x---00000000000000000---x---x---x-
833 : ! 0 | | | 0 | | | 0 |--/| | ! | | H | H | r 0 r | r | r | r 0 e | e | |
834 : ! 0---------------0---------------0--/| |--/| ! x---------------0---------------0---x---x---x-
835 : ! 0 | | | 0 | | | 0 |--/| | ! | | H | H | r 0 r | r | r | r 0 e | e | |
836 : ! 0---------------0---------------0--/| |--/| ! x---------------0---------------0---x---x---x-
837 : ! 0 | | | 0 | | | 0 |--/| 0 ! | | H | H | r 0 r | r | r | r 0 e | e | |
838 : ! 0---------------0---------------0--/| 0000| ! x---------------0---------------0---x---x---x-
839 : ! 0 | | | 0 | | | 0 0000| | ! | | H | H | r 0 r | r | r | r 0 e | e | |
840 : ! 000000000000000000000000000000000000| |--/| ! x---x---x---x---00000000000000000---x---x---x-
841 : ! 0 | | | 0 | | | 0 |--/| | ! | | H | H | r | r | r | r | r | e | e | |
842 : ! 0---------------0---------------0--/| |--/| ! ----------------x---------------x---x---x---x-
843 : ! 0 | | | 0 | | | 0 |--/| | ! | | | H | H | H | H | H | H | e | i | |
844 : ! 0---------------0---------------0--/| |--/| ! ----------------x---------------x---x---x---x-
845 : ! 0 | | | 0 | | | 0 |--/| 0 ! | | | | H | H | H | H | H | i | i | |
846 : ! 0---------------0---------------0--/| 0000 ! ----------------x---------------x---x---x---x-
847 : ! 0 | | | 0 | | | 0 0000 ! | | | | | | | | | i | | |
848 : ! 000000000000000000000000000000000000 ! x---x---x---x---x---x---x---x---x---x---x---x-
849 : !
850 : !
851 : ! -3 |-2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7
852 : !
853 10581322752 : fpanel (1-nht:nc ,1-nht:nc+nht )=fcube(1-nht:nc ,1-nht:nc+nht)
854 9754656912 : w = halo_interp_weight(:,:,:,1)
855 495999504 : do halo=1,nhr
856 2479997520 : do i=halo-nh,nc+nh-(halo-1)
857 1983998016 : ibaseref = ibase(i,halo,1)
858 2314664352 : fpanel (nc+halo ,i ) = matmul_w(w(:,i,halo),fcube(nc +halo,ibaseref:ibaseref+ns-1),ns)
859 : end do
860 : end do
861 :
862 165333168 : if (present(fotherpanel)) then
863 6117327216 : fotherpanel (nc+1 :nc+nht ,1-nht:nc+nht,1)=fcube(nc+1 :nc+nht ,1-nht:nc+nht) !
864 495999504 : do halo=1,nhr
865 2479997520 : do i=halo-nh,nc+nh-(halo-1)
866 : ! ibaseref=ibase(i,halo,1 )
867 1983998016 : ibaseref = ibase(i,halo,1)
868 2314664352 : fotherpanel (nc+1-halo ,i,1) = matmul_w(w(:,i,halo),fcube(nc+1-halo,ibaseref:ibaseref+ns-1),ns)
869 : end do
870 : end do
871 : end if
872 :
873 : else if (cubeboundary==north) then
874 : !
875 : ! north part is on different panel
876 : !
877 : ! stencil
878 : !
879 : ! CN<1 case
880 : ! ! x---------------x---------------x---------------x
881 : ! ! | | | | | | | | | | | | |
882 : !0---\---\---\---0---\---\---\---0---\---\---\---0 ! x---------------x---------------x---------------x
883 : ! 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 ! | | i | i | n | n | n | n | n | n | i | i | |
884 : ! 0---\---\---\---0---\---\---\---0---\---\---\---0 ! x---------------x---------------x---------------x
885 : ! 0 \ \ \ 0 \ \ \ 0 \ \ \ 0 ! | i | i | n | n | n | n | n | n | n | n | i | i |
886 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
887 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
888 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
889 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
890 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
891 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
892 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
893 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
894 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
895 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r | r | r | r | r | r | H | H | |
896 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
897 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | H | H | H | H | H | H | H | H | | |
898 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
899 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | | H | H | H | H | H | H | | | |
900 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
901 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | | | | | | | | | | |
902 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---x---x---x---x---x---x---x---x---x
903 : !
904 : ! -3 |-2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 ! -3 |-2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
905 : !
906 : ! fill in values that are on the same projection as "main" element
907 10085323248 : fpanel (1-nht:nc+nht ,1-nht:nc)=fcube(1-nht:nc+nht ,1-nht:nc)
908 9754656912 : w = halo_interp_weight(:,:,:,1)
909 495999504 : do halo=1,nhr
910 2479997520 : do i=halo-nh,nc+nh-(halo-1)
911 1983998016 : ibaseref = ibase(i,halo,1)
912 2314664352 : fpanel (i,nc+halo ) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+halo ),ns) !north
913 : end do
914 : end do
915 165333168 : if (present(fotherpanel)) then
916 : ! fill in halo for north element
917 5125328208 : fotherpanel (1-nht:nc+nht ,nc+1:nc+nht,1)=fcube(1-nht:nc+nht ,nc+1:nc+nht)
918 : !
919 495999504 : do halo=1,nhr
920 2479997520 : do i=halo-nh,nc+nh-(halo-1)
921 1983998016 : ibaseref = ibase(i,halo,1)
922 2314664352 : fotherpanel (i,nc+1-halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+1-halo),ns)
923 : end do
924 : end do
925 : end if
926 :
927 :
928 : else if (cubeboundary==south) then
929 : !
930 : ! south part is on different panel
931 : !
932 : ! stencil
933 : !
934 : ! !
935 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---x---x---x---x---x---x---x---x---x
936 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | | | | | | | | | | |
937 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
938 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | | H | H | H | H | H | H | | | |
939 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
940 : ! 0 | | | 0 | | | 0 | | | 0 ! | | | H | H | H | H | H | H | H | H | | |
941 : ! 0---------------0---------------0---------------0 ! x---------------x---------------x---------------x
942 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r | r | r | r | r | r | H | H | |
943 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
944 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
945 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
946 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
947 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
948 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
949 : ! 0---------------0---------------0---------------0 ! x---------------0---------------0---------------x
950 : ! 0 | | | 0 | | | 0 | | | 0 ! | | H | H | r 0 r | r | r | r 0 r | H | H | |
951 : ! 0000000000000000000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
952 : ! 0 / / / 0 / / / 0 / / / 0 ! | i | i | s | s | s | s | s | s | s | s | i | i |
953 : ! 0---/---/---/---0---/---/---/---0---/---/---/---0 ! x---------------x---------------x---------------x
954 : ! 0 / / / 0 / / / 0 / / / 0 ! | | i | i | s | s | s | s | s | s | i | i | |
955 : !0---/---/---/---0---/---/---/---0---/---/---/---0 ! x---------------x---------------x---------------x
956 : ! ! | | | | | | | | | | | | |
957 : !
958 : ! 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 ! 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
959 : !
960 : ! fill in values that are on the same projection as "main" element (marked with "i" in Figure above)
961 : !
962 10085323248 : fpanel (1-nht:nc+nht,1:nc+nht )=fcube(1-nht:nc+nht,1:nc+nht)
963 9754656912 : w = halo_interp_weight(:,:,:,1)
964 495999504 : do halo=1,nhr
965 2479997520 : do i=halo-nh,nc+nh-(halo-1)
966 1983998016 : ibaseref=ibase(i,halo,1)!ibase(i,halo,2)
967 2314664352 : fpanel (i,1-halo ) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,1-halo),ns) !south
968 : end do
969 : end do
970 165333168 : if (present(fotherpanel)) then
971 5125328208 : fotherpanel (1-nht:nc+nht,1-nht:0 ,1)=fcube(1-nht:nc+nht,1-nht:0 )
972 495999504 : do halo=1,nhr
973 2479997520 : do i=halo-nh,nc+nh-(halo-1)
974 1983998016 : ibaseref=ibase(i,halo,1)!ibase(i,halo,2)
975 2314664352 : fotherpanel (i, halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1, halo),ns)
976 : end do
977 : end do
978 : end if
979 : else if (cubeboundary==swest) then
980 : !
981 : ! south and west neighboring cells are on different panel
982 : !
983 : ! stencil
984 : !
985 : !
986 : ! CN<1 case
987 : !
988 : !
989 : !
990 : ! |000000000000000000000000000000000000 ! x---x---x---x---x---x---x---x---x---x---x---x---x
991 : ! 0000 0 | | | 0 | | | 0 ! | | | | | | | | | | | | |
992 : ! 0 |/--0---------------0---------------0 ! x---------------x---------------x---------------x
993 : ! |/--| 0 | | | 0 | | | 0 ! | | | | w | H | H | H | H | H | | | |
994 : ! | |/--0---------------0---------------0 ! x---------------x---------------x---------------x
995 : ! |/--| 0 | | | 0 | | | 0 ! | | | w | w | H | H | H | H | H | H | | |
996 : ! | |/--0---------------0---------------0 ! x---------------x---------------x---------------x
997 : ! |/--| 0 | | | 0 | | | 0 ! | | | w | w | r | r | r | r | r | H | H | |
998 : ! | |000000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
999 : ! |0000 0 | | | 0 | | | 0 ! | | | w | w 0 r | r | r | r 0 r | H | H | |
1000 : ! 0 |/--0---------------0---------------0 ! x---------------0---------------0---------------x
1001 : ! |/--| 0 | | | 0 | | | 0 ! | | | w | w 0 r | r | r | r 0 r | H | H | |
1002 : ! | |/--0---------------0---------------0 ! x---------------0---------------0---------------x
1003 : ! |/--| 0 | | | 0 | | | 0 ! | | | w | w 0 r | r | r | r 0 r | H | H | |
1004 : ! | |/--0---------------0---------------0 ! x---------------0---------------0---------------x
1005 : ! | | 0 | | | 0 | | | 0 ! | | | w | w 0 r | r | r | r 0 r | H | H | |
1006 : ! | -/| 000000000000000000000000000000000 ! x---x---x---x---00000000000000000---x---x---x---x
1007 : ! |/ | 0 / / / 0 / / / 0 ! | | | w | | s | s | s | s | s | s | | |
1008 : ! | 0-----/---/---/---0---/---/---/---0 ! x---------------x---------------x---------------x
1009 : ! | 0 / / / 0 / / / 0 ! | | | | s | s | s | s | s | s | | | |
1010 : ! 0-------/---/---/---0---/---/---/---0 ! x---------------x---------------x---------------x
1011 : !
1012 : !
1013 : ! -1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ! |-3 |-2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
1014 : !
1015 : ! fill in values that are on the same projection as "main" element (marked with "i" in Figure above)
1016 : !
1017 253904508 : fpanel(1:nc+nht,1:nc+nht)=fcube(1:nc+nht,1:nc+nht)
1018 : !
1019 : ! fill in west part (marked with "w" on Figure above) and south part (marked with "s")
1020 : !
1021 348380604 : w = halo_interp_weight(:,:,:,1)
1022 17714268 : do halo=1,nhr
1023 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1024 64952316 : ibaseref=ibase(i,halo,1)!ibase(i,halo,1)
1025 64952316 : fpanel(1-halo ,i) = matmul_w(w(:,i,halo),fcube(1-halo ,ibaseref:ibaseref+ns-1),ns) !west
1026 76761828 : fpanel(i,1-halo ) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,1-halo) ,ns) !south
1027 : end do
1028 : end do
1029 : !
1030 : ! corner value
1031 : !
1032 5904756 : fpanel(0,0) =0.25_r8*(fpanel(0,1)+fpanel(1,0)+fpanel(-1,0)+fpanel(0,-1))
1033 : !
1034 : ! ****************************************************************
1035 : !
1036 : ! fill halo for reconstruction on south neighbor panel projection
1037 : !
1038 : ! ****************************************************************
1039 : !
1040 : ! On the south panel projection the neighbors are arragened as follow (nwest case):
1041 : !
1042 : !
1043 : ! \
1044 : ! \ p
1045 : ! \
1046 : ! \-----
1047 : ! |
1048 : ! w | s
1049 : ! |
1050 : !
1051 : !
1052 : ! x---x---x---x---00000000000000000---x---x---x---x
1053 : ! | | | | 0 | | | 0 | | | |
1054 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1055 : ! | | | | 0 | | | 0 | | | |
1056 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1057 : ! | | | | p 0 p | p | p | p 0 p | | | |
1058 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1059 : ! | | | w | wp0 p | p | p | p 0 p | p | | |
1060 : ! x---x---x---x---00000000000000000---x---x---x---x
1061 : ! | | | w | w | r | r | r | r | r | i | i | |
1062 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1063 : ! | | | | w | i | i | i | i | i | i | | |
1064 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1065 : ! | | | | | i | i | i | i | i | | | |
1066 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1067 : ! | | | | | | | | | | | | |
1068 : ! x---x---x---x---x---x---x---x---x---x---x---x---
1069 : !
1070 : !
1071 : ! fill values on same panel projection ("r" and "i" on Figure above)
1072 : !
1073 5904756 : if (present(fotherpanel)) then
1074 129904632 : fotherpanel(1:nc+nht,1-nht:0,1) = fcube(1:nc+nht,1-nht:0)
1075 : !
1076 : ! compute interpolated cell average values in "p" cells on Figure on above
1077 : !
1078 348380604 : w = halo_interp_weight(:,:,:,1)
1079 17714268 : do halo=1,nhr
1080 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1081 64952316 : ibaseref=ibase(i,halo,1)
1082 : !
1083 : ! use same weights as interpolation south from main panel (symmetric)
1084 : !
1085 76761828 : fotherpanel(i,halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,halo),ns)
1086 : end do
1087 : end do
1088 : !
1089 : ! compute interpolated cell average values in "w" cells on Figure on above
1090 : !
1091 348380604 : w = halo_interp_weight(:,:,:,2)
1092 17714268 : do halo=1,nhr
1093 47238048 : do i=nc+halo-nhr,nc+1
1094 29523780 : ibaseref=ibase(i,halo,2)-nc
1095 : !
1096 : ! fotherpanel indexing follows main panel indexing
1097 : ! fcube indexing most be "rotated":
1098 : !
1099 : ! ===============================
1100 : ! | | |
1101 : ! | W ^ | S |
1102 : ! | | | |
1103 : ! | x | | |
1104 : ! | | | |
1105 : ! ! | |
1106 : ! ! <----- | |
1107 : ! ! y | |
1108 : ! ! | |
1109 : ! ===============================
1110 : !
1111 41333292 : fotherpanel(1-halo,i-nc,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,halo),ns)
1112 : end do
1113 : end do
1114 5904756 : fotherpanel(0,1,1) = 0.25_r8*(fotherpanel(-1,1,1)+fotherpanel(1,1,1)+fotherpanel(0,2,1)+fotherpanel(0,0,1))
1115 : !
1116 : ! ****************************************************************
1117 : !
1118 : ! fill halo for reconstruction on west neighbor panel projection
1119 : !
1120 : ! ****************************************************************
1121 : !
1122 : ! On the west panel projection the neighbors are arragened as follow (seast case):
1123 : !
1124 : ! --------
1125 : ! | |
1126 : ! | w | p
1127 : ! | |
1128 : ! -------\
1129 : ! \
1130 : ! s \
1131 : !
1132 : !
1133 : !
1134 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1135 : ! | | | | | | | | | | | | |
1136 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1137 : ! | | | | i | | | | | | | | |
1138 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1139 : ! | | | i | i | e | | | | | | | |
1140 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1141 : ! | | i | i | r | e | e | | | | | | |
1142 : ! x---x---x---x---00000000000000000---x---x---x---x
1143 : ! | | i | i | r 0 e | e | | 0 | | | |
1144 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1145 : ! | | i | i | r 0 e | e | | 0 | | | |
1146 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1147 : ! | | i | i | r 0 e | e | | 0 | | | |
1148 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1149 : ! | | i | i | r 0 e | e | | 0 | | | |
1150 : ! x---x---x---x---00000000000000000---x---x---x---x
1151 : ! | | | s | s | se| e | | | | | | |
1152 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1153 : ! | | | | s | s | | | | | | | |
1154 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1155 : ! | | | | | | | | | | | | |
1156 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1157 : ! | | | | | | | | | | | | |
1158 : ! x---x---x---x---x---x---x---x---x---x---x---x---
1159 : !
1160 : !
1161 : ! fill values on same panel projection ("r" and "i" on Figure above)
1162 : !
1163 253904508 : fotherpanel(1-nht:nc,1:nc+nht,2) = fcube(1-nht:nc,1:nc+nht)
1164 : !
1165 : ! compute interpolated cell average values in "p" cells on Figure on above
1166 : !
1167 348380604 : w = halo_interp_weight(:,:,:,1) ! symmetry
1168 17714268 : do halo=1,nhr
1169 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1170 64952316 : ibaseref=ibase(i,halo,1)
1171 : !
1172 : ! use same weights as interpolation south from main panel (symmetric)
1173 : !
1174 76761828 : fotherpanel(halo,i,2) = matmul_w(w(:,i,halo),fcube(halo,ibaseref:ibaseref+ns-1),ns)
1175 : end do
1176 : end do
1177 : !
1178 : ! compute interpolated cell average values in "s" cells on Figure on above
1179 : !
1180 348380604 : w = halo_interp_weight(:,:,:,2)
1181 17714268 : do halo=1,nhr
1182 47238048 : do i=nc+halo-nhr,nc+1
1183 29523780 : ibaseref=ibase(i,halo,2)-nc
1184 : !
1185 : ! fotherpanel indexing follows main panel indexing
1186 : ! fcube indexing most be "rotated":
1187 : !
1188 : ! ===============================
1189 : ! | | |
1190 : ! | W ^ | S |
1191 : ! | | | |
1192 : ! | x | | |
1193 : ! | | | |
1194 : ! ! | |
1195 : ! ! <----- | |
1196 : ! ! y | |
1197 : ! ! | |
1198 : ! ===============================
1199 : !
1200 41333292 : fotherpanel(i-nc,1-halo,2) = matmul_w(w(:,i,halo),fcube(halo,ibaseref:ibaseref+ns-1),ns)
1201 : end do
1202 : end do
1203 5904756 : fotherpanel(1,0,2) = 0.25_r8*(fotherpanel(0,0,2)+fotherpanel(2,0,2)+fotherpanel(1,-1,2)+fotherpanel(1,1,2))
1204 : end if
1205 : else if (cubeboundary==seast) then
1206 : !
1207 : ! south and east neighboring cells are on different panel
1208 : !
1209 : !
1210 : !
1211 : ! 000000000000000000000000000000000000|
1212 : ! 0 | | | 0 | | | 0 0000 ! | | | | | | | | | | | | |
1213 : ! 0---------------0---------------0--\| 0 ! x---------------x---------------x---------------x
1214 : ! 0 | | | 0 | | | 0 |--\| ! | | | | | H | H | H | H | | | | |
1215 : ! 0---------------0---------------0--\| | ! x---------------x---------------x---------------x
1216 : ! 0 | | | 0 | | | 0 |--\| ! | | | | H | H | H | H | H | e | | | |
1217 : ! 0---------------0---------------0--\| | ! x---------------x---------------x---------------x
1218 : ! 0 | | | 0 | | | 0 |--\| ! | | H | H | r | r | r | r | r | e | e | | |
1219 : ! 000000000000000000000000000000000000| | ! x---x---x---x---00000000000000000---x---x---x---x
1220 : ! 0 | | | 0 | | | 0 0000| ! | | H | H | r 0 r | r | r | r 0 e | e | | |
1221 : ! 0---------------0---------------0--\| 0 ! x---------------0---------------0---------------x
1222 : ! 0 | | | 0 | | | 0 |--\| ! | | H | H | r 0 r | r | r | r 0 e | e | | |
1223 : ! 0---------------0---------------0--\| | ! x---------------0---------------0---------------x
1224 : ! 0 | | | 0 | | | 0 |--\| ! | | H | H | r 0 r | r | r | r 0 e | e | | |
1225 : ! 0---------------0---------------0--\| | ! x---------------0---------------0---------------x
1226 : ! 0 | | | 0 | | | 0 | | ! | | H | H | r 0 r | r | r | r 0 e | e | | |
1227 : ! 000000000000000000000000000000000 |\- | ! x---x---x---x---00000000000000000---x---x---x---x
1228 : ! 0 \ \ \ 0 \ \ \ 0 | \| ! | | | s | s | s | s | s | s |s/e| e | | |
1229 : ! 0---\---\---\---0---\---\---\-----0 | ! x---------------x---------------x---------------x
1230 : ! 0 \ \ \ 0 \ \ \ 0 | ! | | | | s | s | s | s | s | s | | | |
1231 : ! 0---\---\---\---0---\---\---\-------0 ! x---------------x---------------x---------------x
1232 : !
1233 : !
1234 253904508 : fpanel (1-nht:nc,1:nc+nht)=fcube(1-nht:nc,1:nc+nht)
1235 : !
1236 : ! east
1237 : !
1238 348380604 : w = halo_interp_weight(:,:,:,1)
1239 17714268 : do halo=1,nhr
1240 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1241 64952316 : ibaseref = ibase(i,halo,1)
1242 76761828 : fpanel(nc+halo,i) = matmul_w(w(:,i,halo),fcube(nc +halo,ibaseref:ibaseref+ns-1),ns)
1243 : end do
1244 : end do
1245 : !
1246 : ! south
1247 : !
1248 348380604 : w = halo_interp_weight(:,:,:,2)
1249 17714268 : do halo=1,nhr
1250 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1251 64952316 : ibaseref = ibase(i,halo,2)
1252 76761828 : fpanel(i,1-halo ) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,1-halo),ns) !south
1253 : end do
1254 : end do
1255 11809512 : fpanel(nc+1,0 )=0.25_r8*(&
1256 11809512 : fpanel(nc+1,1)+fpanel(nc,0)+fpanel(nc+2,0)+fpanel(nc+1,-1))
1257 : !
1258 : ! ****************************************************************
1259 : !
1260 : ! fill halo for reconstruction on south neighbor panel projection
1261 : !
1262 : ! ****************************************************************
1263 : !
1264 : ! On the south panel projection the neighbors are arragened as follow (neast case):
1265 : !
1266 : !
1267 : ! /
1268 : ! P /
1269 : ! /
1270 : ! ------/
1271 : ! | | E
1272 : ! | S |
1273 : ! | |
1274 : !
1275 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1276 : ! | | | | | | | | | | | | |
1277 : ! x---x---x---x---00000000000000000---x---x---x---x
1278 : ! | | | | 0 | | | 0 | | | |
1279 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1280 : ! | | | | 0 | | | 0 | | | |
1281 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1282 : ! | | | | n 0 n | n | n | n 0 n | | | |
1283 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1284 : ! | | | n | n 0 n | n | n | n 0 ne| e | | |
1285 : ! x---x---x---x---00000000000000000---x---x---x---x
1286 : ! | | i | i | r | r | r | r | r | e | e | | |
1287 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1288 : ! | | | i | i | i | i | i | i | e | | | |
1289 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1290 : ! | | | | i | i | i | i | i | | | | |
1291 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1292 : ! | | | | | | | | | | | | |
1293 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1294 : !
1295 : !
1296 : !
1297 : ! fill values on same panel projection ("r" and "i" on Figure above)
1298 : !
1299 5904756 : if (present(fotherpanel)) then
1300 129904632 : fotherpanel(1-nht:nc,1-nht:0,1) = fcube(1-nht:nc,1-nht:0)
1301 : !
1302 348380604 : w = halo_interp_weight(:,:,:,2)
1303 : !
1304 : ! fill in "n" on Figure above
1305 : !
1306 17714268 : do halo=1,nhr
1307 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1308 64952316 : ibaseref = ibase(i,halo,2)
1309 76761828 : fotherpanel (i,halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1, halo),ns)
1310 : end do
1311 : end do
1312 : !
1313 : ! fill in "e" on Figure above
1314 : !
1315 348380604 : w = halo_interp_weight(:,:,:,1)
1316 17714268 : do halo=1,nhr
1317 47238048 : do i=0,nht-halo!nc+nh-(halo-1)
1318 29523780 : ibaseref = ibase(i,halo,1)
1319 : !
1320 : ! fother panel follows indexing on main panel
1321 : !
1322 : ! use symmetry for weights (same weights as East from main panel but for south panel
1323 : ! projection the indecies are rotated)
1324 : !
1325 41333292 : fotherpanel (nc+halo ,1-i,1) = matmul_w(w(:,i,halo),fcube(nc+ibaseref:nc+ibaseref+ns-1,halo),ns)
1326 : end do
1327 : end do
1328 5904756 : fotherpanel(nc+1,1,1) = 0.25_r8*(fotherpanel(nc+2,1,1)+fotherpanel(nc,1,1)&
1329 5904756 : +fotherpanel(nc+1,2,1)+fotherpanel(nc+1,0,1))
1330 :
1331 : !
1332 : ! ****************************************************************
1333 : !
1334 : ! fill halo for reconstruction on east neighbor panel projection
1335 : !
1336 : ! ****************************************************************
1337 : !
1338 : ! On the south panel projection the neighbors are arragened as follow (neast case):
1339 : !
1340 : !
1341 : ! | |
1342 : ! P | E |
1343 : ! |-----|
1344 : ! /
1345 : ! / S
1346 : ! /
1347 : !
1348 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1349 : ! | | | | | | | | | | | | |
1350 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1351 : ! | | | | | | | | | i | | | |
1352 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1353 : ! | | | | | | | | w | i | i | | |
1354 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1355 : ! | | | | | | | w | w | r | i | i | |
1356 : ! x---x---x---x---00000000000000000---x---x---x---x
1357 : ! | | | | 0 | | w | w 0 r | i | i | |
1358 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1359 : ! | | | | 0 | | w | w 0 r | i | i | |
1360 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1361 : ! | | | | 0 | | w | w 0 r | i | i | |
1362 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1363 : ! | | | | 0 | | w | w 0 r | i | i | |
1364 : ! x---x---x---x---00000000000000000---x---x---x---x
1365 : ! | | | | | | | w | ws| s | s | | |
1366 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1367 : ! | | | | | | | | s | s | | | |
1368 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1369 : ! | | | | | | | | | | | | |
1370 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1371 : !
1372 : !
1373 : !
1374 : ! fill values on same panel projection ("r" and "i" on Figure above)
1375 : !
1376 147618900 : fotherpanel(nc+1:nc+nht,1:nc+nht,2) = fcube(nc+1:nc+nht,1:nc+nht)
1377 : !
1378 : !
1379 : ! fill in "w" on Figure above
1380 : !
1381 348380604 : w = halo_interp_weight(:,:,:,1)
1382 17714268 : do halo=1,nhr
1383 82666584 : do i=0,nc+nh-(halo-1)
1384 64952316 : ibaseref = ibase(i,halo,1)
1385 76761828 : fotherpanel(nc+1-halo,i,2) = matmul_w(w(:,i,halo),fcube(nc+1-halo,ibaseref:ibaseref+ns-1),ns)
1386 : end do
1387 : end do
1388 : !
1389 : ! fill in "s" on Figure above
1390 : !
1391 348380604 : w = halo_interp_weight(:,:,:,2)
1392 17714268 : do halo=1,nhr
1393 47238048 : do i=nc+1-nht+halo,nc+1
1394 : !
1395 : !
1396 : ! ! P | E
1397 : ! ! |
1398 : ! ! |
1399 : ! ================
1400 : ! | |
1401 : ! | S | | <----- y
1402 : ! | | | ^
1403 : ! | x | | |
1404 : ! | v | |
1405 : ! ! | |
1406 : ! ! -----> | x
1407 : ! ! y |
1408 : ! ! |
1409 : ! ================
1410 : !
1411 : !
1412 : ! shift (since we are using south weights from main panel interpolation
1413 : !
1414 29523780 : ibaseref = ibase(i,halo,2)-nc
1415 : !
1416 : ! fotherpanel index: reverse
1417 : !
1418 : ! fcube index: due to rotation (see Figure above)
1419 : !
1420 41333292 : fotherpanel(nc+(nc+1-i),1-halo,2) = matmul_w(w(:,i,halo),fcube(nc+1-halo,ibaseref:ibaseref+ns-1),ns)
1421 : end do
1422 : end do
1423 5904756 : fotherpanel(nc,0,2) = 0.25_r8*(fotherpanel(nc+1,0,2)+fotherpanel(nc-1,0,2)&
1424 11809512 : +fotherpanel(nc,1,2)+fotherpanel(nc,-1,2))
1425 : end if
1426 : else if (cubeboundary==nwest) then
1427 : !
1428 : !
1429 : ! 0-------\---\---\---0---\---\---\---0 ! --------x---------------x---------------x
1430 : ! | 0 \ \ \ 0 \ \ \ 0 ! | | n | n | n | n | n | n | | | |
1431 : ! | 0-----\---\---\---0---\---\---\---0 ! --------x---------------x---------------x
1432 : ! | | 0 \ \ \ 0 \ \ \ 0 ! | w | a | n | n | n | n | n | n | | |
1433 : ! |\ | 000000000000000000000000000000000 ! --------00000000000000000---------------x
1434 : ! | -\| 0 | | | 0 | | | 0 ! | w | w 0 r | r | r | r 0 r | H | H | |
1435 : ! | |\--0---------------0---------------0 ! --------0---------------0---------------x
1436 : ! |\--| 0 | | | 0 | | | 0 ! | w | w 0 r | r | r | r 0 r | H | H | |
1437 : ! | |\--0---------------0---------------0 ! --------0---------------0---------------x
1438 : ! |\--| 0 | | | 0 | | | 0 ! | w | w 0 r | r | r | r 0 r | H | H | |
1439 : ! 0 |\--0---------------0---------------0 ! --------0---------------0---------------x
1440 : ! |0000 0 | | | 0 | | | 0 ! | w | w 0 r | r | r | r 0 r | H | H | |
1441 : ! | |000000000000000000000000000000000000 ! --------00000000000000000---------------x
1442 : ! |\--| 0 | | | 0 | | | 0 ! | w | w | r | r | r | r | r | H | H | |
1443 : ! | |\--0---------------0---------------0 ! --------x---------------x---------------x
1444 : ! |\--| 0 | | | 0 | | | 0 ! | | w | H | H | H | H | H | H | | |
1445 : ! | |\--0---------------0---------------0 ! --------x---------------x---------------x
1446 : ! |\--| 0 | | | 0 | | | 0 ! | | | H | H | H | H | H | | | |
1447 : ! 0 |\--0---------------0---------------0 ! --------x---------------x---------------x
1448 : ! 0000 0 | | | 0 | | | 0 ! | | | | | | | | | | |
1449 : ! 000000000000000000000000000000000000 ! --------x---------------x---------------x
1450 : !
1451 : !
1452 : !
1453 253904508 : fpanel(1:nc+nht,1-nht:nc)=fcube(1:nc+nht,1-nht:nc)
1454 : !
1455 : ! west
1456 : !
1457 348380604 : w = halo_interp_weight(:,:,:,1)
1458 17714268 : do halo=1,nhr
1459 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1460 64952316 : ibaseref=ibase(i,halo,1)
1461 76761828 : fpanel(1-halo ,i) = matmul_w(w(:,i,halo),fcube(1-halo ,ibaseref:ibaseref+ns-1),ns)
1462 : end do
1463 : end do
1464 : !
1465 : ! north
1466 : !
1467 348380604 : w = halo_interp_weight(:,:,:,2)
1468 17714268 : do halo=1,nhr
1469 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1470 64952316 : ibaseref = ibase(i,halo,2)
1471 76761828 : fpanel(i,nc+halo) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+halo ),ns) !north
1472 : end do
1473 : end do
1474 5904756 : fpanel(0 ,nc+1)=0.25_r8*(&
1475 11809512 : fpanel(0,nc)+fpanel(1,nc+1)+fpanel(-1,nc+1)+fpanel(0,nc+2))
1476 : !
1477 : ! ****************************************************************
1478 : !
1479 : ! fill halo for reconstruction on north neighbor panel projection
1480 : !
1481 : ! ****************************************************************
1482 : !
1483 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1484 : !| | | | | | | | | | | | |
1485 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1486 : !| | | | | i | i | i | i | i | | | |
1487 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1488 : !| | | | w | i | i | i | i | i | i | | |
1489 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1490 : !| | | w | w | r | r | r | r | r | i | i | |
1491 : !x---x---x---x---00000000000000000---x---x---x---x
1492 : !| | | w | ws0 s | s | s | s 0 s | s | | |
1493 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1494 : !| | | | s 0 s | s | s | s 0 s | | | |
1495 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1496 : !| | | | 0 | | | 0 | | | |
1497 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1498 : !| | | | 0 | | | 0 | | | |
1499 : !x---x---x---x---00000000000000000---x---x---x---x
1500 : !
1501 : !
1502 : ! fill values on same panel projection ("r" and "i" on Figure above)
1503 : !
1504 5904756 : if (present(fotherpanel)) then
1505 129904632 : fotherpanel(1:nc+nht,nc+1:nc+nht,1) = fcube(1:nc+nht,nc+1:nc+nht)
1506 : !
1507 : !
1508 : ! fill in "s" on Figure above
1509 : !
1510 : ! (use code from north above)
1511 : !
1512 348380604 : w = halo_interp_weight(:,:,:,2)
1513 17714268 : do halo=1,nhr
1514 82666584 : do i=max(halo-nh,0),nc+nh-(halo-1)
1515 64952316 : ibaseref = ibase(i,halo,2)
1516 76761828 : fotherpanel(i,nc+1-halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+1-halo ),ns)
1517 : end do
1518 : end do
1519 : !
1520 : ! fill in "w" on Figure above
1521 : !
1522 : ! (use code from west above)
1523 : !
1524 348380604 : w = halo_interp_weight(:,:,:,1)
1525 17714268 : do halo=1,nhr
1526 47238048 : do i=nc+1-nht+halo,nc+1
1527 29523780 : ibaseref=ibase(i,halo,1)-nc
1528 41333292 : fotherpanel(1-halo,nc-(i-(nc+1)),1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+1-halo),ns)
1529 : end do
1530 : end do
1531 5904756 : fotherpanel(0,nc,1)=0.25_r8*(&
1532 5904756 : fotherpanel(1,nc,1)+fotherpanel(-1,nc,1)+fotherpanel(0,nc+1,1)+fotherpanel(0,nc-1,1))
1533 :
1534 : !
1535 : ! ****************************************************************
1536 : !
1537 : ! fill halo for reconstruction on west neighbor panel projection
1538 : !
1539 : ! ****************************************************************
1540 : !
1541 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1542 : !| | | | | | | | | | | | |
1543 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1544 : !| | | | | | | | | | | | |
1545 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1546 : !| | | | n | n | | | | | | | |
1547 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1548 : !| | | n | n | ne| e | | | | | | |
1549 : !x---x---x---x---00000000000000000---x---x---x---x
1550 : !| | i | i | r 0 e | e | | 0 | | | |
1551 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1552 : !| | i | i | r 0 e | e | | 0 | | | |
1553 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1554 : !| | i | i | r 0 e | e | | 0 | | | |
1555 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1556 : !| | i | i | r 0 e | e | | 0 | | | |
1557 : !x---x---x---x---00000000000000000---x---x---x---x
1558 : !| | i | i | r | e | e | | | | | | |
1559 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1560 : !| | | i | i | e | | | | | | | |
1561 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1562 : !| | | | i | | | | | | | | |
1563 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1564 : !| | | | | | | | | | | | |
1565 : !x---x---x---x---x---x---x---x---x---x---x---x---
1566 : !
1567 : !
1568 : ! fill values on same panel projection ("r" and "i" on Figure above)
1569 : !
1570 253904508 : fotherpanel(1-nht:nc,1-nht:nc,2) = fcube(1-nht:nc,1-nht:nc)
1571 : !
1572 : !
1573 : ! fill in "e" on Figure above
1574 : !
1575 : ! (use code from west above)
1576 : !
1577 348380604 : w = halo_interp_weight(:,:,:,1)
1578 17714268 : do halo=1,nhr
1579 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1580 64952316 : ibaseref=ibase(i,halo,1)
1581 76761828 : fotherpanel(halo ,i,2) = matmul_w(w(:,i,halo),fcube(halo ,ibaseref:ibaseref+ns-1),ns)
1582 : end do
1583 : end do
1584 : !
1585 : !
1586 : ! fill in "n" on Figure above
1587 : !
1588 : ! (use code from north above)
1589 : !
1590 348380604 : w = halo_interp_weight(:,:,:,2)
1591 17714268 : do halo=1,nhr
1592 47238048 : do i=0,nht-halo
1593 29523780 : ibaseref = ibase(i,halo,2)+nc
1594 41333292 : fotherpanel(1-i,nc+halo,2) = matmul_w(w(:,i,halo),fcube(halo,ibaseref:ibaseref+ns-1),ns) !north
1595 : end do
1596 : end do
1597 5904756 : fotherpanel(1,nc+1,2)=0.25_r8*(&
1598 5904756 : fotherpanel(2,nc+1,2)+fotherpanel(0,nc+1,2)+fotherpanel(1,nc+2,2)+fotherpanel(1,nc,2))
1599 : end if
1600 :
1601 : else if (cubeboundary==neast) then
1602 : !
1603 : !
1604 : ! 0---/---/---/---0---/---/---/-------0 ! x---------------x---------------x--------
1605 : ! 0 / / / 0 / / / 0 | ! | | | | | n | n | n | n | n | |
1606 : ! 0---/---/---/---0---/---/---/-----0 | ! x---------------x---------------x--------
1607 : ! 0 / / / 0 / / / 0 | | ! | | | | n | n | n | n | n | a | e |
1608 : ! 000000000000000000000000000000000 | | ! x---------------00000000000000000--------
1609 : ! 0 | | | 0 | | | 0 |--/| ! | | | H | H 0 r | r | r | r 0 e | e |
1610 : ! 0---------------0---------------0--/| | ! x---------------0---------------0--------
1611 : ! 0 | | | 0 | | | 0 |--/| ! | | | H | H 0 r | r | r | r 0 e | e |
1612 : ! 0---------------0---------------0--/| | ! x---------------0---------------0--------
1613 : ! 0 | | | 0 | | | 0 |--/| ! | | | H | H 0 r | r | r | r 0 e | e |
1614 : ! 0---------------0---------------0--/| 0 ! x---------------0---------------0--------
1615 : ! 0 | | | 0 | | | 0 0000| ! | | | H | H 0 r | r | r | r 0 e | e |
1616 : ! 000000000000000000000000000000000000| | ! x---------------00000000000000000--------
1617 : ! 0 | | | 0 | | | 0 |--/| ! | | | H | H | r | r | r | r | e | e |
1618 : ! 0---------------0---------------0--/| | ! x---------------x---------------x--------
1619 : ! 0 | | | 0 | | | 0 |--/| ! | | | | H | H | H | H | H | e | |
1620 : ! 0---------------0---------------0--/| | ! x---------------x---------------x--------
1621 : ! 0 | | | 0 | | | 0 |--/| ! | | | | | H | H | H | H | | |
1622 : ! 0---------------0---------------0--/| 0 ! x---------------x---------------x--------
1623 : ! 0 | | | 0 | | | 0 0000 ! | | | | | | | | | | |
1624 : ! 000000000000000000000000000000000000 ! x---------------x---------------x--------
1625 : !
1626 : !
1627 : !
1628 253904508 : fpanel(1-nht:nc,1-nht:nc)=fcube(1-nht:nc,1-nht:nc)
1629 : ! fotherpanel (nc+1 :nc+nht ,1-nht:nc+nht)=fcube(nc+1 :nc+nht ,1-nht:nc+nht)
1630 : !
1631 : ! east
1632 : !
1633 348380604 : w = halo_interp_weight(:,:,:,1)
1634 17714268 : do halo=1,nhr
1635 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1636 64952316 : ibaseref=ibase(i,halo,1 )
1637 76761828 : fpanel(nc+halo,i) = matmul_w(w(:,i,halo),fcube(nc +halo,ibaseref:ibaseref+ns-1),ns)
1638 : end do
1639 : end do
1640 : !
1641 : ! north
1642 : !
1643 : ! w = halo_interp_weight(:,:,:,1)
1644 17714268 : do halo=1,nhr
1645 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1646 64952316 : ibaseref=ibase(i,halo,1)
1647 76761828 : fpanel(i,nc+halo) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+halo ),ns) !north
1648 : end do
1649 : end do
1650 5904756 : fpanel(nc+1,nc+1)=0.25_r8*(&
1651 11809512 : fpanel(nc,nc+1)+fpanel(nc+1,nc)+fpanel(nc+1,nc+2)+fpanel(nc+2,nc+1))
1652 : !
1653 : ! ****************************************************************
1654 : !
1655 : ! fill halo for reconstruction on north neighbor panel projection
1656 : !
1657 : ! ****************************************************************
1658 : !
1659 : ! On the north panel projection the neighbors are arragened as follow (seast case):
1660 : !
1661 : !
1662 : ! | |
1663 : ! | N | E
1664 : ! |-----|
1665 : ! \
1666 : ! S \
1667 : ! \
1668 : !
1669 : !
1670 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1671 : ! | | | | | | | | | | | | |
1672 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1673 : ! | | | | i | i | i | i | i | | | | |
1674 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1675 : ! | | | i | i | i | i | i | i | e | | | |
1676 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1677 : ! | | i | i | r | r | r | r | r | e | e | | |
1678 : ! x---x---x---x---00000000000000000---x---x---x---x
1679 : ! | | | s | s 0 s | s | s | s 0 se| e | | |
1680 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1681 : ! | | | | s 0 s | s | s | s 0 s | | | |
1682 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1683 : ! | | | | 0 | | | 0 | | | |
1684 : ! x---x---x---x---0---x---x---x---0---x---x---x---x
1685 : ! | | | | 0 | | | 0 | | | |
1686 : ! x---x---x---x---00000000000000000---x---x---x---x
1687 : ! | | | | | | | | | | | | |
1688 : ! x---x---x---x---x---x---x---x---x---x---x---x---x
1689 : !
1690 : !
1691 : ! fill values on same panel projection ("r" and "i" on Figure above)
1692 : !
1693 5904756 : if (present(fotherpanel)) then
1694 129904632 : fotherpanel(1-nht:nc,nc+1:nc+nht,1) = fcube(1-nht:nc,nc+1:nc+nht)
1695 : !
1696 : ! fill in "s" on Figure above
1697 : !
1698 : ! (use north case from above and shift/reverse j-index
1699 : !
1700 348380604 : w = halo_interp_weight(:,:,:,1)
1701 17714268 : do halo=1,nhr
1702 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1703 64952316 : ibaseref=ibase(i,halo,1)
1704 76761828 : fotherpanel (i,nc+1-halo,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+1-halo),ns)
1705 : end do
1706 : end do
1707 : !
1708 : ! fill in "e" on Figure above
1709 : !
1710 348380604 : w = halo_interp_weight(:,:,:,2)
1711 17714268 : do halo=1,nhr
1712 47238048 : do i=max(halo-nh,0),nht-halo
1713 29523780 : ibaseref=ibase(i,halo,2) +nc
1714 : !
1715 : ! fotherpanel uses indexing of main panel's projection
1716 : ! fcube: rotated indexing
1717 : !
1718 41333292 : fotherpanel (nc+halo,nc+i,1) = matmul_w(w(:,i,halo),fcube(ibaseref:ibaseref+ns-1,nc+1-halo),ns)
1719 : end do
1720 : end do
1721 5904756 : fotherpanel(nc+1,nc,1)=0.25_r8*(&
1722 5904756 : fotherpanel(nc+2,nc,1)+fotherpanel(nc,nc,1)+fotherpanel(nc+1,nc+1,1)+fotherpanel(nc+1,nc-1,1))
1723 : !
1724 : ! ****************************************************************
1725 : !
1726 : ! fill halo for reconstruction on east neighbor panel projection
1727 : !
1728 : ! ****************************************************************
1729 : !
1730 : ! On the north panel projection the neighbors are arragened as follow (seast case):
1731 : !
1732 : !
1733 : ! \ N
1734 : ! \
1735 : ! \------
1736 : ! | |
1737 : ! P | E |
1738 : ! | |
1739 : ! -------
1740 : !
1741 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1742 : !| | | | | | | | | | | | |
1743 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1744 : !| | | | | | | | | | | | |
1745 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1746 : !| | | | | | | | n | n | | | |
1747 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1748 : !| | | | | | | w | wn| n | n | | |
1749 : !x---x---x---x---00000000000000000---x---x---x---x
1750 : !| | | | 0 | | w | w 0 r | i | i | |
1751 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1752 : !| | | | 0 | | w | w 0 r | i | i | |
1753 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1754 : !| | | | 0 | | w | w 0 r | i | i | |
1755 : !x---x---x---x---0---x---x---x---0---x---x---x---x
1756 : !| | | | 0 | | w | w 0 r | i | i | |
1757 : !x---x---x---x---00000000000000000---x---x---x---x
1758 : !| | | | | | | w | w | r | i | i | |
1759 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1760 : !| | | | | | | | w | i | i | | |
1761 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1762 : !| | | | | | | | | i | | | |
1763 : !x---x---x---x---x---x---x---x---x---x---x---x---x
1764 : !| | | | | | | | | | | | |
1765 : !x---x---x---x---x---x---x---x---x---x---x---x---
1766 : !
1767 : !
1768 : !
1769 : ! fill values on same panel projection ("r" and "i" on Figure above)
1770 : !
1771 147618900 : fotherpanel(nc+1:nc+nht,1-nht:nc,2) = fcube(nc+1:nc+nht,1-nht:nc)
1772 : !
1773 : ! fill in "w" on Figure above
1774 : !
1775 : ! (use east case from above and shift/reverse j-index
1776 : !
1777 348380604 : w = halo_interp_weight(:,:,:,1)
1778 17714268 : do halo=1,nhr
1779 82666584 : do i=halo-nh,min(nc+nh-(halo-1),nc+1)
1780 64952316 : ibaseref=ibase(i,halo,1 )
1781 76761828 : fotherpanel(nc+1-halo,i,2) = matmul_w(w(:,i,halo),fcube(nc+1-halo,ibaseref:ibaseref+ns-1),ns)
1782 : end do
1783 : end do
1784 : !
1785 : ! fill in "n" on Figure above
1786 : !
1787 348380604 : w = halo_interp_weight(:,:,:,2)
1788 17714268 : do halo=1,nhr
1789 47238048 : do i=max(halo-nh,0),nht-halo
1790 29523780 : ibaseref=ibase(i,halo,2) +nc
1791 : !
1792 : ! fotherpanel uses indexing of main panel's projection
1793 : ! fcube: rotated indexing
1794 : !
1795 41333292 : fotherpanel (nc+i,nc+halo,2) = matmul_w(w(:,i,halo),fcube(nc+1-halo,ibaseref:ibaseref+ns-1),ns)
1796 : end do
1797 : end do
1798 : fotherpanel(nc,nc+1,2)=0.25_r8*(&
1799 5904756 : fotherpanel(nc+1,nc+1,2)+fotherpanel(nc-1,nc+1,2)+fotherpanel(nc,nc+2,2)+fotherpanel(nc,nc,2))
1800 : end if
1801 : end if
1802 5314280400 : end subroutine extend_panel_interpolate
1803 : !
1804 : ! initialize non-existent ghost cells
1805 : !
1806 483116400 : subroutine zero_non_existent_ghost_cell(recons,irecons,cubeboundary,nc,nhe,ntrac_in)
1807 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
1808 :
1809 : integer, intent(in) :: nc,nhe,cubeboundary,irecons,ntrac_in
1810 : real (kind=r8), dimension(irecons,1-nhe:nc+nhe,1-nhe:nc+nhe,ntrac_in), intent(out):: recons
1811 :
1812 : integer :: i,j
1813 :
1814 483116400 : if (cubeboundary>0) then
1815 62268336 : if (cubeboundary==nwest) then
1816 1073592 : do j=nc+1,nc+nhe
1817 1610388 : do i=1-nhe,0
1818 42406884 : recons(:,i,j,:) = 0.0_r8
1819 : end do
1820 : end do
1821 61731540 : else if (cubeboundary==swest) then
1822 1073592 : do j=1-nhe,0
1823 1610388 : do i=1-nhe,0
1824 42406884 : recons(:,i,j,:) = 0.0_r8
1825 : end do
1826 : end do
1827 61194744 : else if (cubeboundary==neast) then
1828 1073592 : do j=nc+1,nc+nhe
1829 1610388 : do i=nc+1,nc+nhe
1830 42406884 : recons(:,i,j,:) = 0.0_r8
1831 : end do
1832 : end do
1833 60657948 : else if (cubeboundary==seast) then
1834 1073592 : do j=1-nhe,0
1835 1610388 : do i=nc+1,nc+nhe
1836 42406884 : recons(:,i,j,:) = 0.0_r8
1837 : end do
1838 : end do
1839 : end if
1840 : end if
1841 483116400 : end subroutine zero_non_existent_ghost_cell
1842 : end module fvm_reconstruction_mod
|
