Line data Source code
1 : module cube_mod
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 : use coordinate_systems_mod, only: spherical_polar_t, cartesian3D_t, cartesian2d_t, &
4 : projectpoint, cubedsphere2cart, spherical_to_cart, sphere_tri_area,dist_threshold, &
5 : change_coordinates
6 :
7 : use physconst, only: pi, rearth
8 : use control_mod, only: hypervis_scaling, cubed_sphere_map
9 : use cam_abortutils, only: endrun
10 :
11 : implicit none
12 : private
13 :
14 : integer,public, parameter :: nfaces = 6 ! number of faces on the cube
15 : integer,public, parameter :: nInnerElemEdge = 8 ! number of edges for an interior element
16 : integer,public, parameter :: nCornerElemEdge = 4 ! number of corner elements
17 :
18 : real(kind=r8), public, parameter :: cube_xstart = -0.25_R8*PI
19 : real(kind=r8), public, parameter :: cube_xend = 0.25_R8*PI
20 : real(kind=r8), public, parameter :: cube_ystart = -0.25_R8*PI
21 : real(kind=r8), public, parameter :: cube_yend = 0.25_R8*PI
22 :
23 :
24 : type, public :: face_t
25 : type (spherical_polar_t) :: sphere0 ! tangent point of face on sphere
26 : type (spherical_polar_t) :: sw ! sw corner of face on sphere
27 : type (spherical_polar_t) :: se ! se corner of face on sphere
28 : type (spherical_polar_t) :: ne ! ne corner of face on sphere
29 : type (spherical_polar_t) :: nw ! nw corner of face on sphere
30 : type (cartesian3D_t) :: P0
31 : type (cartesian3D_t) :: X0
32 : type (cartesian3D_t) :: Y0
33 : integer :: number
34 : integer :: padding ! pad the struct
35 : end type face_t
36 :
37 : type, public :: cube_face_coord_t
38 : real(r8) :: x ! x coordinate
39 : real(r8) :: y ! y coordinate
40 : type (face_t), pointer :: face ! face
41 : end type cube_face_coord_t
42 :
43 : ! ==========================================
44 : ! Public Interfaces
45 : ! ==========================================
46 :
47 : public :: CubeTopology
48 :
49 : ! Rotate the North Pole: used for JW baroclinic test case
50 : ! Settings this only changes Coriolis.
51 : ! User must also rotate initial condition
52 : real (kind=r8), public :: rotate_grid = 0
53 :
54 : ! ===============================
55 : ! Public methods for cube
56 : ! ===============================
57 :
58 : public :: cube_init_atomic
59 : public :: convert_gbl_index
60 : public :: vmap,dmap
61 : public :: covariant_rot
62 : public :: contravariant_rot
63 : public :: set_corner_coordinates
64 : public :: assign_node_numbers_to_elem
65 :
66 :
67 : public :: CubeEdgeCount
68 : public :: CubeElemCount
69 : public :: CubeSetupEdgeIndex
70 : public :: rotation_init_atomic
71 : public :: ref2sphere
72 :
73 : ! ===============================
74 : ! Private methods
75 : ! ===============================
76 : private :: coordinates_atomic
77 : private :: metric_atomic
78 : private :: coreolis_init_atomic
79 :
80 : contains
81 :
82 : ! =======================================
83 : ! cube_init_atomic:
84 : !
85 : ! Initialize element descriptors for
86 : ! cube sphere case for each element ...
87 : ! =======================================
88 21600 : subroutine cube_init_atomic(elem, gll_points, alpha_in)
89 : use element_mod, only : element_t
90 : use dimensions_mod, only : np
91 :
92 : type (element_t), intent(inout) :: elem
93 : real(r8), intent(in) :: gll_points(np)
94 : real(r8), optional, intent(in) :: alpha_in
95 :
96 : real(r8) :: alpha
97 :
98 21600 : if (present(alpha_in)) then
99 10800 : alpha = alpha_in
100 : else
101 10800 : alpha = 1.0_r8
102 : end if
103 :
104 21600 : elem%FaceNum=elem%vertex%face_number
105 21600 : call coordinates_atomic(elem,gll_points)
106 :
107 21600 : call metric_atomic(elem, gll_points, alpha)
108 :
109 21600 : call coreolis_init_atomic(elem)
110 21600 : elem%desc%use_rotation= 0
111 :
112 21600 : end subroutine cube_init_atomic
113 :
114 : ! =======================================
115 : ! coordinates_atomic:
116 : !
117 : ! Initialize element coordinates for
118 : ! cube-sphere case ... (atomic)
119 : !
120 : ! =======================================
121 :
122 21600 : subroutine coordinates_atomic(elem, gll_points)
123 : use element_mod, only: element_t, element_var_coordinates
124 : use dimensions_mod, only: np
125 :
126 : type(element_t), intent(inout) :: elem
127 : real(r8), intent(in) :: gll_points(np)
128 :
129 : real(r8) :: area1,area2
130 : type (cartesian3d_t) :: quad(4)
131 : integer :: face_no,i,j
132 :
133 21600 : face_no = elem%vertex%face_number
134 : ! compute the corners in Cartesian coordinates
135 108000 : do i=1,4
136 108000 : elem%corners3D(i)=cubedsphere2cart(elem%corners(i),face_no)
137 : enddo
138 :
139 : ! =========================================
140 : ! compute lat/lon coordinates of each GLL point
141 : ! =========================================
142 108000 : do i=1,np
143 453600 : do j=1,np
144 432000 : elem%spherep(i,j)=ref2sphere(gll_points(i),gll_points(j),elem%corners3D,cubed_sphere_map,elem%corners,elem%facenum)
145 : enddo
146 : enddo
147 :
148 : ! also compute the [-pi/2,pi/2] cubed sphere coordinates:
149 453600 : elem%cartp=element_var_coordinates(elem%corners,gll_points)
150 :
151 : ! Matrix describing vector conversion to cartesian
152 : ! Zonal direction
153 453600 : elem%vec_sphere2cart(:,:,1,1) = -SIN(elem%spherep(:,:)%lon)
154 453600 : elem%vec_sphere2cart(:,:,2,1) = COS(elem%spherep(:,:)%lon)
155 453600 : elem%vec_sphere2cart(:,:,3,1) = 0.0_r8
156 : ! Meridional direction
157 453600 : elem%vec_sphere2cart(:,:,1,2) = -SIN(elem%spherep(:,:)%lat)*COS(elem%spherep(:,:)%lon)
158 453600 : elem%vec_sphere2cart(:,:,2,2) = -SIN(elem%spherep(:,:)%lat)*SIN(elem%spherep(:,:)%lon)
159 453600 : elem%vec_sphere2cart(:,:,3,2) = COS(elem%spherep(:,:)%lat)
160 :
161 21600 : end subroutine coordinates_atomic
162 :
163 : ! elem_jacobians:
164 : !
165 : ! Calculate Jacobian associated with mapping
166 : ! from arbitrary quadrilateral to [-1,1]^2
167 : ! along with its inverse and determinant
168 : ! ==========================================
169 :
170 21600 : subroutine elem_jacobians(coords, unif2quadmap)
171 :
172 : use dimensions_mod, only : np
173 : type (cartesian2D_t), dimension(np,np), intent(in) :: coords
174 : ! unif2quadmap is the bilinear map from [-1,1]^2 -> arbitrary quadrilateral
175 : real (kind=r8), dimension(4,2), intent(out) :: unif2quadmap
176 : integer :: ii,jj
177 :
178 21600 : unif2quadmap(1,1)=(coords(1,1)%x+coords(np,1)%x+coords(np,np)%x+coords(1,np)%x)/4.0_r8
179 21600 : unif2quadmap(1,2)=(coords(1,1)%y+coords(np,1)%y+coords(np,np)%y+coords(1,np)%y)/4.0_r8
180 21600 : unif2quadmap(2,1)=(-coords(1,1)%x+coords(np,1)%x+coords(np,np)%x-coords(1,np)%x)/4.0_r8
181 21600 : unif2quadmap(2,2)=(-coords(1,1)%y+coords(np,1)%y+coords(np,np)%y-coords(1,np)%y)/4.0_r8
182 21600 : unif2quadmap(3,1)=(-coords(1,1)%x-coords(np,1)%x+coords(np,np)%x+coords(1,np)%x)/4.0_r8
183 21600 : unif2quadmap(3,2)=(-coords(1,1)%y-coords(np,1)%y+coords(np,np)%y+coords(1,np)%y)/4.0_r8
184 21600 : unif2quadmap(4,1)=(coords(1,1)%x-coords(np,1)%x+coords(np,np)%x-coords(1,np)%x)/4.0_r8
185 21600 : unif2quadmap(4,2)=(coords(1,1)%y-coords(np,1)%y+coords(np,np)%y-coords(1,np)%y)/4.0_r8
186 :
187 21600 : end subroutine elem_jacobians
188 :
189 : ! =========================================
190 : ! metric_atomic:
191 : !
192 : ! Initialize cube-sphere metric terms:
193 : ! equal angular elements (atomic)
194 : ! initialize:
195 : ! metdet, rmetdet (analytic) = detD, 1/detD
196 : ! met (analytic) D^t D (symmetric)
197 : ! metdet (analytic) = detD
198 : ! metinv (analytic) Dinv Dinv^t (symmetic)
199 : ! D (from subroutine vmap)
200 : ! Dinv (computed directly from D)
201 : !
202 : ! ucontra = Dinv * u = metinv * ucov
203 : ! ucov = D^t * u = met * ucontra
204 : !
205 : ! we also compute DE = D*E, where
206 : ! E = eigenvectors of metinv as a basis metinv = E LAMBDA E^t
207 : !
208 : ! ueig = E^t ucov = E^t D^t u = (DE)^t u
209 : !
210 : !
211 : ! so if we want to tweak the mapping by a factor alpha (so he weights add up to 4pi, for example)
212 : ! we take:
213 : ! NEW OLD
214 : ! D = sqrt(alpha) D and then rederive all quantities.
215 : ! detD = alpha detD
216 : !
217 : ! where alpha = 4pi/SEMarea, SEMarea = global sum elem(ie)%mv(i,j)*elem(ie)%metdet(i,j)
218 : !
219 : ! =========================================
220 :
221 21600 : subroutine metric_atomic(elem,gll_points,alpha)
222 : use element_mod, only: element_t
223 : use dimensions_mod, only: np
224 : use physconst, only: ra
225 :
226 : type (element_t), intent(inout) :: elem
227 : real(r8), intent(in) :: alpha
228 : real(r8), intent(in) :: gll_points(np)
229 :
230 : ! Local variables
231 : integer ii
232 : integer i,j,nn
233 : integer iptr
234 :
235 : real (kind=r8) :: r ! distance from origin for point on cube tangent to unit sphere
236 :
237 : real (kind=r8) :: const, norm
238 : real (kind=r8) :: detD ! determinant of vector field mapping matrix.
239 :
240 : real (kind=r8) :: x1 ! 1st cube face coordinate
241 : real (kind=r8) :: x2 ! 2nd cube face coordinate
242 : real (kind=r8) :: tmpD(2,2)
243 : real (kind=r8) :: M(2,2),E(2,2),eig(2),DE(2,2),DEL(2,2),V(2,2), nu1, nu2, lamStar1, lamStar2
244 : integer :: imaxM(2)
245 : real (kind=r8) :: l1, l2, sc,min_svd,max_svd,max_normDinv
246 :
247 :
248 : ! ==============================================
249 : ! Initialize differential mapping operator
250 : ! to and from vector fields on the sphere to
251 : ! contravariant vector fields on the cube
252 : ! i.e. dM/dx^i in Sadourney (1972) and it's
253 : ! inverse
254 : ! ==============================================
255 :
256 : ! MNL: Calculate Jacobians of bilinear map from cubed-sphere to ref element
257 21600 : if (cubed_sphere_map==0) then
258 21600 : call elem_jacobians(elem%cartp, elem%u2qmap)
259 : endif
260 :
261 : max_svd = 0.0_r8
262 : max_normDinv = 0.0_r8
263 : min_svd = 1d99
264 108000 : do j=1,np
265 453600 : do i=1,np
266 345600 : x1=gll_points(i)
267 345600 : x2=gll_points(j)
268 2419200 : call Dmap(elem%D(i,j,:,:),x1,x2,elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
269 :
270 :
271 : ! Numerical metric tensor based on analytic D: met = D^T times D
272 : ! (D maps between sphere and reference element)
273 : elem%met(i,j,1,1) = elem%D(i,j,1,1)*elem%D(i,j,1,1) + &
274 345600 : elem%D(i,j,2,1)*elem%D(i,j,2,1)
275 : elem%met(i,j,1,2) = elem%D(i,j,1,1)*elem%D(i,j,1,2) + &
276 345600 : elem%D(i,j,2,1)*elem%D(i,j,2,2)
277 : elem%met(i,j,2,1) = elem%D(i,j,1,1)*elem%D(i,j,1,2) + &
278 345600 : elem%D(i,j,2,1)*elem%D(i,j,2,2)
279 : elem%met(i,j,2,2) = elem%D(i,j,1,2)*elem%D(i,j,1,2) + &
280 345600 : elem%D(i,j,2,2)*elem%D(i,j,2,2)
281 :
282 : ! compute D^-1...
283 : ! compute determinant of D mapping matrix... if not zero compute inverse
284 :
285 345600 : detD = elem%D(i,j,1,1)*elem%D(i,j,2,2) - elem%D(i,j,1,2)*elem%D(i,j,2,1)
286 :
287 345600 : elem%Dinv(i,j,1,1) = elem%D(i,j,2,2)/detD
288 345600 : elem%Dinv(i,j,1,2) = -elem%D(i,j,1,2)/detD
289 345600 : elem%Dinv(i,j,2,1) = -elem%D(i,j,2,1)/detD
290 345600 : elem%Dinv(i,j,2,2) = elem%D(i,j,1,1)/detD
291 :
292 : ! L2 norm = sqrt max eigenvalue of metinv
293 : ! = 1/sqrt(min eigenvalue of met)
294 : ! l1 and l2 are eigenvalues of met
295 : ! (should both be positive, l1 > l2)
296 : l1 = (elem%met(i,j,1,1) + elem%met(i,j,2,2) + sqrt(4.0_r8*elem%met(i,j,1,2)*elem%met(i,j,2,1) + &
297 345600 : (elem%met(i,j,1,1) - elem%met(i,j,2,2))**2))/2.0_r8
298 : l2 = (elem%met(i,j,1,1) + elem%met(i,j,2,2) - sqrt(4.0_r8*elem%met(i,j,1,2)*elem%met(i,j,2,1) + &
299 345600 : (elem%met(i,j,1,1) - elem%met(i,j,2,2))**2))/2.0_r8
300 : ! Max L2 norm of Dinv is sqrt of max eigenvalue of metinv
301 : ! max eigenvalue of metinv is 1/min eigenvalue of met
302 345600 : norm = 1.0_r8/sqrt(min(abs(l1),abs(l2)))
303 345600 : max_svd = max(norm, max_svd)
304 : ! Min L2 norm of Dinv is sqrt of min eigenvalue of metinv
305 : ! min eigenvalue of metinv is 1/max eigenvalue of met
306 345600 : norm = 1.0_r8/sqrt(max(abs(l1),abs(l2)))
307 345600 : min_svd = min(norm, min_svd)
308 :
309 : ! some kind of pseudo-norm of Dinv
310 : ! C = 1/sqrt(2) sqrt( |g^x|^2 + |g^y|^2 + 2*|g^x dot g^y|)
311 : ! = 1/sqrt(2) sqrt( |g_x|^2 + |g_y|^2 + 2*|g_x dot g_y|) / J
312 : ! g^x = Dinv(:,1) g_x = D(1,:)
313 : ! g^y = Dinv(:,2) g_y = D(2,:)
314 2419200 : norm = (2*abs(sum(elem%Dinv(i,j,:,1)*elem%Dinv(i,j,:,2))) + sum(elem%Dinv(i,j,:,1)**2) + sum(elem%Dinv(i,j,:,2)**2))
315 345600 : norm = sqrt(norm)
316 : ! norm = (2*abs(sum(elem%D(1,:,i,j)*elem%D(2,:,i,j))) + sum(elem%D(1,:,i,j)**2) + sum(elem%D(2,:,i,j)**2))
317 : ! norm = sqrt(norm)/detD
318 : max_normDinv = max(norm,max_normDinv)
319 :
320 :
321 : ! Need inverse of met if not calculated analytically
322 345600 : elem%metdet(i,j) = abs(detD)
323 345600 : elem%rmetdet(i,j) = 1.0_R8/abs(detD)
324 :
325 345600 : elem%metinv(i,j,1,1) = elem%met(i,j,2,2)/(detD*detD)
326 345600 : elem%metinv(i,j,1,2) = -elem%met(i,j,1,2)/(detD*detD)
327 345600 : elem%metinv(i,j,2,1) = -elem%met(i,j,2,1)/(detD*detD)
328 345600 : elem%metinv(i,j,2,2) = elem%met(i,j,1,1)/(detD*detD)
329 :
330 : ! matricies for tensor hyper-viscosity
331 : ! compute eigenvectors of metinv (probably same as computed above)
332 2419200 : M = elem%metinv(i,j,:,:)
333 :
334 : eig(1) = (M(1,1) + M(2,2) + sqrt(4.0_r8*M(1,2)*M(2,1) + &
335 345600 : (M(1,1) - M(2,2))**2))/2.0_r8
336 : eig(2) = (M(1,1) + M(2,2) - sqrt(4.0_r8*M(1,2)*M(2,1) + &
337 345600 : (M(1,1) - M(2,2))**2))/2.0_r8
338 :
339 : ! use DE to store M - Lambda, to compute eigenvectors
340 345600 : DE=M
341 345600 : DE(1,1)=DE(1,1)-eig(1)
342 345600 : DE(2,2)=DE(2,2)-eig(1)
343 :
344 2419200 : imaxM = maxloc(abs(DE))
345 2419200 : if (maxval(abs(DE))==0) then
346 48 : E(1,1)=1; E(2,1)=0;
347 345552 : elseif ( imaxM(1)==1 .and. imaxM(2)==1 ) then
348 172784 : E(2,1)=1; E(1,1) = -DE(2,1)/DE(1,1)
349 172768 : else if ( imaxM(1)==1 .and. imaxM(2)==2 ) then
350 0 : E(2,1)=1; E(1,1) = -DE(2,2)/DE(1,2)
351 172768 : else if ( imaxM(1)==2 .and. imaxM(2)==1 ) then
352 912 : E(1,1)=1; E(2,1) = -DE(1,1)/DE(2,1)
353 171856 : else if ( imaxM(1)==2 .and. imaxM(2)==2 ) then
354 171856 : E(1,1)=1; E(2,1) = -DE(1,2)/DE(2,2)
355 : else
356 0 : call endrun('Impossible error in cube_mod.F90::metric_atomic()')
357 : endif
358 :
359 : ! the other eigenvector is orthgonal:
360 345600 : E(1,2)=-E(2,1)
361 345600 : E(2,2)= E(1,1)
362 :
363 : !normalize columns
364 1728000 : E(:,1)=E(:,1)/sqrt(sum(E(:,1)*E(:,1)));
365 1728000 : E(:,2)=E(:,2)/sqrt(sum(E(:,2)*E(:,2)));
366 :
367 :
368 : ! OBTAINING TENSOR FOR HV:
369 :
370 : ! Instead of the traditional scalar Laplace operator \grad \cdot \grad
371 : ! we introduce \grad \cdot V \grad
372 : ! where V = D E LAM LAM^* E^T D^T.
373 : ! Recall (metric_tensor)^{-1}=(D^T D)^{-1} = E LAM E^T.
374 : ! Here, LAM = diag( 4/((np-1)dx)^2 , 4/((np-1)dy)^2 ) = diag( 4/(dx_elem)^2, 4/(dy_elem)^2 )
375 : ! Note that metric tensors and LAM correspondingly are quantities on a unit sphere.
376 :
377 : ! This motivates us to use V = D E LAM LAM^* E^T D^T
378 : ! where LAM^* = diag( nu1, nu2 ) where nu1, nu2 are HV coefficients scaled like (dx)^{hv_scaling/2}, (dy)^{hv_scaling/2}.
379 : ! (Halves in powers come from the fact that HV consists of two Laplace iterations.)
380 :
381 : ! Originally, we took LAM^* = diag(
382 : ! 1/(eig(1)**(hypervis_scaling/4.0_r8))*(rearth**(hypervis_scaling/2.0_r8))
383 : ! 1/(eig(2)**(hypervis_scaling/4.0_r8))*(rearth**(hypervis_scaling/2.0_r8)) ) =
384 : ! = diag( lamStar1, lamStar2)
385 : ! \simeq ((np-1)*dx_sphere / 2 )^hv_scaling/2 = SQRT(OPERATOR_HV)
386 : ! because 1/eig(...) \simeq (dx_on_unit_sphere)^2 .
387 : ! Introducing the notation OPERATOR = lamStar^2 is useful for conversion formulas.
388 :
389 : ! This leads to the following conversion formula: nu_const is nu used for traditional HV on uniform grids
390 : ! nu_tensor = nu_const * OPERATOR_HV^{-1}, so
391 : ! nu_tensor = nu_const *((np-1)*dx_sphere / 2 )^{ - hv_scaling} or
392 : ! nu_tensor = nu_const *(2/( (np-1) * dx_sphere) )^{hv_scaling} .
393 : ! dx_sphere = 2\pi *rearth/(np-1)/4/NE
394 : ! [nu_tensor] = [meter]^{4-hp_scaling}/[sec]
395 :
396 : ! (1) Later developments:
397 : ! Apply tensor V only at the second Laplace iteration. Thus, LAM^* should be scaled as (dx)^{hv_scaling}, (dy)^{hv_scaling},
398 : ! see this code below:
399 : ! DEL(1:2,1) = (lamStar1**2) *eig(1)*DE(1:2,1)
400 : ! DEL(1:2,2) = (lamStar2**2) *eig(2)*DE(1:2,2)
401 :
402 : ! (2) Later developments:
403 : ! Bringing [nu_tensor] to 1/[sec]:
404 : ! lamStar1=1/(eig(1)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
405 : ! lamStar2=1/(eig(2)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
406 : ! OPERATOR_HV = ( (np-1)*dx_unif_sphere / 2 )^{hv_scaling} * rearth^4
407 : ! Conversion formula:
408 : ! nu_tensor = nu_const * OPERATOR_HV^{-1}, so
409 : ! nu_tensor = nu_const *( 2*rearth /((np-1)*dx))^{hv_scaling} * rearth^{-4.0}.
410 :
411 : ! For the baseline coefficient nu=1e15 for NE30,
412 : ! nu_tensor=7e-8 (BUT RUN TWICE AS SMALL VALUE FOR NOW) for hv_scaling=3.2
413 : ! and
414 : ! nu_tensor=1.3e-6 for hv_scaling=4.0.
415 :
416 :
417 : !matrix D*E
418 1036800 : DE(1,1)=sum(elem%D(i,j,1,:)*E(:,1))
419 1036800 : DE(1,2)=sum(elem%D(i,j,1,:)*E(:,2))
420 1036800 : DE(2,1)=sum(elem%D(i,j,2,:)*E(:,1))
421 1036800 : DE(2,2)=sum(elem%D(i,j,2,:)*E(:,2))
422 :
423 345600 : lamStar1=1/(eig(1)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
424 345600 : lamStar2=1/(eig(2)**(hypervis_scaling/4.0_r8)) *(rearth**2.0_r8)
425 :
426 : !matrix (DE) * Lam^* * Lam , tensor HV when V is applied at each Laplace calculation
427 : ! DEL(1:2,1) = lamStar1*eig(1)*DE(1:2,1)
428 : ! DEL(1:2,2) = lamStar2*eig(2)*DE(1:2,2)
429 :
430 : !matrix (DE) * (Lam^*)^2 * Lam, tensor HV when V is applied only once, at the last Laplace calculation
431 : !will only work with hyperviscosity, not viscosity
432 1036800 : DEL(1:2,1) = (lamStar1**2) *eig(1)*DE(1:2,1)
433 1036800 : DEL(1:2,2) = (lamStar2**2) *eig(2)*DE(1:2,2)
434 :
435 : !matrix (DE) * Lam^* * Lam *E^t *D^t or (DE) * (Lam^*)^2 * Lam *E^t *D^t
436 1036800 : V(1,1)=sum(DEL(1,:)*DE(1,:))
437 1036800 : V(1,2)=sum(DEL(1,:)*DE(2,:))
438 1036800 : V(2,1)=sum(DEL(2,:)*DE(1,:))
439 1036800 : V(2,2)=sum(DEL(2,:)*DE(2,:))
440 :
441 2505600 : elem%tensorVisc(i,j,:,:)=V(:,:)
442 :
443 : end do
444 : end do
445 :
446 : ! see Paul Ullrich writeup:
447 : ! max_normDinv might be a tighter bound than max_svd for deformed elements
448 : ! max_svd >= max_normDinv/sqrt(2), with equality holding if |g^x| = |g^y|
449 : ! elem%normDinv=max_normDinv/sqrt(2)
450 :
451 : ! this norm is consistent with length scales defined below:
452 21600 : elem%normDinv=max_svd
453 :
454 :
455 : ! compute element length scales, based on SVDs, in km:
456 21600 : elem%dx_short = 1.0_r8/(max_svd*0.5_r8*dble(np-1)*ra*1000.0_r8)
457 21600 : elem%dx_long = 1.0_r8/(min_svd*0.5_r8*dble(np-1)*ra*1000.0_r8)
458 :
459 : ! optional noramlization:
460 1879200 : elem%D = elem%D * sqrt(alpha)
461 1879200 : elem%Dinv = elem%Dinv / sqrt(alpha)
462 453600 : elem%metdet = elem%metdet * alpha
463 453600 : elem%rmetdet = elem%rmetdet / alpha
464 1879200 : elem%met = elem%met * alpha
465 1879200 : elem%metinv = elem%metinv / alpha
466 :
467 21600 : end subroutine metric_atomic
468 :
469 :
470 : ! ========================================
471 : ! covariant_rot:
472 : !
473 : ! 2 x 2 matrix multiply: Db^T * Da^-T
474 : ! for edge rotations: maps face a to face b
475 : !
476 : ! ========================================
477 :
478 0 : function covariant_rot(Da,Db) result(R)
479 :
480 : real (kind=r8) :: Da(2,2)
481 : real (kind=r8) :: Db(2,2)
482 : real (kind=r8) :: R(2,2)
483 :
484 : real (kind=r8) :: detDa
485 :
486 0 : detDa = Da(2,2)*Da(1,1) - Da(1,2)*Da(2,1)
487 :
488 0 : R(1,1)=(Da(2,2)*Db(1,1) - Da(1,2)*Db(2,1))/detDa
489 0 : R(1,2)=(Da(1,1)*Db(2,1) - Da(2,1)*Db(1,1))/detDa
490 0 : R(2,1)=(Da(2,2)*Db(1,2) - Da(1,2)*Db(2,2))/detDa
491 0 : R(2,2)=(Da(1,1)*Db(2,2) - Da(2,1)*Db(1,2))/detDa
492 :
493 0 : end function covariant_rot
494 :
495 : ! ========================================
496 : ! contravariant_rot:
497 : !
498 : ! 2 x 2 matrix multiply: Db^-1 * Da
499 : ! that maps a contravariant vector field
500 : ! from an edge of cube face a to a contiguous
501 : ! edge of cube face b.
502 : !
503 : ! ========================================
504 :
505 8544 : function contravariant_rot(Da,Db) result(R)
506 :
507 : real(kind=r8), intent(in) :: Da(2,2)
508 : real(kind=r8), intent(in) :: Db(2,2)
509 : real(kind=r8) :: R(2,2)
510 :
511 : real(kind=r8) :: detDb
512 :
513 8544 : detDb = Db(2,2)*Db(1,1) - Db(1,2)*Db(2,1)
514 :
515 8544 : R(1,1)=(Da(1,1)*Db(2,2) - Da(2,1)*Db(1,2))/detDb
516 8544 : R(1,2)=(Da(1,2)*Db(2,2) - Da(2,2)*Db(1,2))/detDb
517 8544 : R(2,1)=(Da(2,1)*Db(1,1) - Da(1,1)*Db(2,1))/detDb
518 8544 : R(2,2)=(Da(2,2)*Db(1,1) - Da(1,2)*Db(2,1))/detDb
519 :
520 8544 : end function contravariant_rot
521 :
522 : ! ========================================================
523 : ! Dmap:
524 : !
525 : ! Initialize mapping that tranforms contravariant
526 : ! vector fields on the reference element onto vector fields on
527 : ! the sphere.
528 : ! ========================================================
529 24694200 : subroutine Dmap(D, a,b, corners3D, ref_map, corners, u2qmap, facenum)
530 : real (kind=r8), intent(out) :: D(2,2)
531 : real (kind=r8), intent(in) :: a,b
532 : type (cartesian3D_t) :: corners3D(4) !x,y,z coords of element corners
533 : integer :: ref_map
534 : ! only needed for ref_map=0,1
535 : type (cartesian2D_t),optional :: corners(4) ! gnomonic coords of element corners
536 : real (kind=r8),optional :: u2qmap(4,2)
537 : integer,optional :: facenum
538 :
539 :
540 :
541 24694200 : if (ref_map==0) then
542 24694200 : if (.not. present ( corners ) ) &
543 0 : call endrun('Dmap(): missing arguments for equiangular map')
544 24694200 : call dmap_equiangular(D,a,b,corners,u2qmap,facenum)
545 0 : else if (ref_map==1) then
546 0 : call endrun('equi-distance gnomonic map not yet implemented')
547 0 : else if (ref_map==2) then
548 0 : call dmap_elementlocal(D,a,b,corners3D)
549 : else
550 0 : call endrun('bad value of ref_map')
551 : endif
552 24694200 : end subroutine Dmap
553 :
554 :
555 :
556 : ! ========================================================
557 : ! Dmap:
558 : !
559 : ! Equiangular Gnomonic Projection
560 : ! Composition of equiangular Gnomonic projection to cubed-sphere face,
561 : ! followd by bilinear map to reference element
562 : ! ========================================================
563 24694200 : subroutine dmap_equiangular(D, a,b, corners,u2qmap,facenum )
564 : use dimensions_mod, only : np
565 : real (kind=r8), intent(out) :: D(2,2)
566 : real (kind=r8), intent(in) :: a,b
567 : real (kind=r8) :: u2qmap(4,2)
568 : type (cartesian2D_t) :: corners(4) ! gnomonic coords of element corners
569 : integer :: facenum
570 : ! local
571 : real (kind=r8) :: tmpD(2,2), Jp(2,2),x1,x2,pi,pj,qi,qj
572 : real (kind=r8), dimension(4,2) :: unif2quadmap
573 :
574 : #if 0
575 : ! we shoud get rid of elem%u2qmap() and routine cube_mod.F90::elem_jacobian()
576 : ! and replace with this code below:
577 : ! but this produces roundoff level changes
578 : !unif2quadmap(1,1)=(elem%cartp(1,1)%x+elem%cartp(np,1)%x+elem%cartp(np,np)%x+elem%cartp(1,np)%x)/4.0_r8
579 : !unif2quadmap(1,2)=(elem%cartp(1,1)%y+elem%cartp(np,1)%y+elem%cartp(np,np)%y+elem%cartp(1,np)%y)/4.0_r8
580 : unif2quadmap(2,1)=(-elem%cartp(1,1)%x+elem%cartp(np,1)%x+elem%cartp(np,np)%x-elem%cartp(1,np)%x)/4.0_r8
581 : unif2quadmap(2,2)=(-elem%cartp(1,1)%y+elem%cartp(np,1)%y+elem%cartp(np,np)%y-elem%cartp(1,np)%y)/4.0_r8
582 : unif2quadmap(3,1)=(-elem%cartp(1,1)%x-elem%cartp(np,1)%x+elem%cartp(np,np)%x+elem%cartp(1,np)%x)/4.0_r8
583 : unif2quadmap(3,2)=(-elem%cartp(1,1)%y-elem%cartp(np,1)%y+elem%cartp(np,np)%y+elem%cartp(1,np)%y)/4.0_r8
584 : unif2quadmap(4,1)=(elem%cartp(1,1)%x-elem%cartp(np,1)%x+elem%cartp(np,np)%x-elem%cartp(1,np)%x)/4.0_r8
585 : unif2quadmap(4,2)=(elem%cartp(1,1)%y-elem%cartp(np,1)%y+elem%cartp(np,np)%y-elem%cartp(1,np)%y)/4.0_r8
586 : Jp(1,1) = unif2quadmap(2,1) + unif2quadmap(4,1)*b
587 : Jp(1,2) = unif2quadmap(3,1) + unif2quadmap(4,1)*a
588 : Jp(2,1) = unif2quadmap(2,2) + unif2quadmap(4,2)*b
589 : Jp(2,2) = unif2quadmap(3,2) + unif2quadmap(4,2)*a
590 : #else
591 : ! input (a,b) shold be a point in the reference element [-1,1]
592 : ! compute Jp(a,b)
593 24694200 : Jp(1,1) = u2qmap(2,1) + u2qmap(4,1)*b
594 24694200 : Jp(1,2) = u2qmap(3,1) + u2qmap(4,1)*a
595 24694200 : Jp(2,1) = u2qmap(2,2) + u2qmap(4,2)*b
596 24694200 : Jp(2,2) = u2qmap(3,2) + u2qmap(4,2)*a
597 : #endif
598 :
599 : ! map (a,b) to the [-pi/2,pi/2] equi angular cube face: x1,x2
600 : ! a = gp%points(i)
601 : ! b = gp%points(j)
602 24694200 : pi = (1-a)/2
603 24694200 : pj = (1-b)/2
604 24694200 : qi = (1+a)/2
605 24694200 : qj = (1+b)/2
606 : x1 = pi*pj*corners(1)%x &
607 : + qi*pj*corners(2)%x &
608 : + qi*qj*corners(3)%x &
609 24694200 : + pi*qj*corners(4)%x
610 : x2 = pi*pj*corners(1)%y &
611 : + qi*pj*corners(2)%y &
612 : + qi*qj*corners(3)%y &
613 24694200 : + pi*qj*corners(4)%y
614 :
615 24694200 : call vmap(tmpD,x1,x2,facenum)
616 :
617 : ! Include map from element -> ref element in D
618 24694200 : D(1,1) = tmpD(1,1)*Jp(1,1) + tmpD(1,2)*Jp(2,1)
619 24694200 : D(1,2) = tmpD(1,1)*Jp(1,2) + tmpD(1,2)*Jp(2,2)
620 24694200 : D(2,1) = tmpD(2,1)*Jp(1,1) + tmpD(2,2)*Jp(2,1)
621 24694200 : D(2,2) = tmpD(2,1)*Jp(1,2) + tmpD(2,2)*Jp(2,2)
622 24694200 : end subroutine dmap_equiangular
623 :
624 :
625 :
626 : ! ========================================================
627 : ! vmap:
628 : !
629 : ! Initialize mapping that tranforms contravariant
630 : ! vector fields on the cube onto vector fields on
631 : ! the sphere. This follows Taylor's D matrix
632 : !
633 : ! | cos(theta)dlambda/dx1 cos(theta)dlambda/dx2 |
634 : ! D = | |
635 : ! | dtheta/dx1 dtheta/dx2 |
636 : !
637 : ! ========================================================
638 :
639 24711288 : subroutine vmap(D, x1, x2, face_no)
640 : real(kind=r8), intent(inout) :: D(2,2)
641 : real(kind=r8), intent(in) :: x1
642 : real(kind=r8), intent(in) :: x2
643 : integer, intent(in) :: face_no
644 :
645 : ! Local variables
646 :
647 : real (kind=r8) :: poledist ! SQRT(TAN(x1)**2 +TAN(x2)**2)
648 : real (kind=r8) :: r ! distance from cube point to center of sphere
649 :
650 : real (kind=r8) :: D11
651 : real (kind=r8) :: D12
652 : real (kind=r8) :: D21
653 : real (kind=r8) :: D22
654 : character(len=64) :: errmsg
655 :
656 24711288 : r = SQRT(1.0_r8 + TAN(x1)**2 + TAN(x2)**2)
657 :
658 24711288 : if (face_no >= 1 .and. face_no <= 4) then
659 :
660 16474192 : D11 = 1.0_r8 / (r * COS(x1))
661 16474192 : D12 = 0.0_r8
662 16474192 : D21 = -TAN(x1)*TAN(x2) / (COS(x1)*r*r)
663 16474192 : D22 = 1.0_r8 / (r*r*COS(x1)*COS(x2)*COS(x2))
664 :
665 16474192 : D(1,1) = D11
666 16474192 : D(1,2) = D12
667 16474192 : D(2,1) = D21
668 16474192 : D(2,2) = D22
669 :
670 :
671 8237096 : else if (face_no == 6) then
672 4118548 : poledist = SQRT( TAN(x1)**2 + TAN(x2)**2)
673 4118548 : if (poledist <= DIST_THRESHOLD) then
674 :
675 : ! we set the D transform to the identity matrix
676 : ! which works ONLY for swtc1, phi starting at
677 : ! 3*PI/2... assumes lon at pole == 0
678 :
679 16 : D(1,1) = 1.0_r8
680 16 : D(1,2) = 0.0_r8
681 16 : D(2,1) = 0.0_r8
682 16 : D(2,2) = 1.0_r8
683 :
684 : else
685 :
686 4118532 : D11 = -TAN(x2)/(poledist*COS(x1)*COS(x1)*r)
687 4118532 : D12 = TAN(x1)/(poledist*COS(x2)*COS(x2)*r)
688 4118532 : D21 = -TAN(x1)/(poledist*COS(x1)*COS(x1)*r*r)
689 4118532 : D22 = -TAN(x2)/(poledist*COS(x2)*COS(x2)*r*r)
690 :
691 4118532 : D(1,1) = D11
692 4118532 : D(1,2) = D12
693 4118532 : D(2,1) = D21
694 4118532 : D(2,2) = D22
695 :
696 : end if
697 4118548 : else if (face_no == 5) then
698 4118548 : poledist = SQRT( TAN(x1)**2 + TAN(x2)**2)
699 4118548 : if (poledist <= DIST_THRESHOLD) then
700 :
701 : ! we set the D transform to the identity matrix
702 : ! which works ONLY for swtc1, phi starting at
703 : ! 3*PI/2... assumes lon at pole == 0, i.e. very specific
704 :
705 16 : D(1,1) = 1.0_r8
706 16 : D(1,2) = 0.0_r8
707 16 : D(2,1) = 0.0_r8
708 16 : D(2,2) = 1.0_r8
709 :
710 : else
711 :
712 4118532 : D11 = TAN(x2)/(poledist*COS(x1)*COS(x1)*r)
713 4118532 : D12 = -TAN(x1)/(poledist*COS(x2)*COS(x2)*r)
714 4118532 : D21 = TAN(x1)/(poledist*COS(x1)*COS(x1)*r*r)
715 4118532 : D22 = TAN(x2)/(poledist*COS(x2)*COS(x2)*r*r)
716 :
717 4118532 : D(1,1) = D11
718 4118532 : D(1,2) = D12
719 4118532 : D(2,1) = D21
720 4118532 : D(2,2) = D22
721 :
722 : end if
723 : else
724 0 : write(errmsg, '(a,i0)') 'VMAP: Bad face number, ',face_no
725 0 : call endrun(errmsg)
726 : end if
727 :
728 24711288 : end subroutine vmap
729 :
730 :
731 :
732 :
733 : ! ========================================================
734 : ! Dmap:
735 : !
736 : ! Initialize mapping that tranforms contravariant
737 : ! vector fields on the reference element onto vector fields on
738 : ! the sphere.
739 : ! For Gnomonic, followed by bilinear, this code uses the old vmap()
740 : ! for unstructured grids, this code uses the parametric map that
741 : ! maps quads on the sphere directly to the reference element
742 : ! ========================================================
743 0 : subroutine dmap_elementlocal(D, a,b, corners3D)
744 : use element_mod, only : element_t
745 :
746 : type (element_t) :: elem
747 : real (kind=r8), intent(out) :: D(2,2)
748 : real (kind=r8), intent(in) :: a,b
749 : type (cartesian3d_t) :: corners3D(4)
750 :
751 : type (spherical_polar_t) :: sphere
752 :
753 : real(kind=r8) :: c(3,4), q(4), xx(3), r, lam, th, dd(4,2)
754 : real(kind=r8) :: sinlam, sinth, coslam, costh
755 : real(kind=r8) :: D1(2,3), D2(3,3), D3(3,2), D4(3,2)
756 : integer :: i,j
757 :
758 0 : sphere = ref2sphere(a,b,corners3D,2) ! use element local map, ref_map=2
759 :
760 0 : c(1,1)=corners3D(1)%x; c(2,1)=corners3D(1)%y; c(3,1)=corners3D(1)%z;
761 0 : c(1,2)=corners3D(2)%x; c(2,2)=corners3D(2)%y; c(3,2)=corners3D(2)%z;
762 0 : c(1,3)=corners3D(3)%x; c(2,3)=corners3D(3)%y; c(3,3)=corners3D(3)%z;
763 0 : c(1,4)=corners3D(4)%x; c(2,4)=corners3D(4)%y; c(3,4)=corners3D(4)%z;
764 :
765 0 : q(1)=(1-a)*(1-b); q(2)=(1+a)*(1-b); q(3)=(1+a)*(1+b); q(4)=(1-a)*(1+b);
766 0 : q=q/4.0_r8;
767 :
768 0 : do i=1,3
769 0 : xx(i)=sum(c(i,:)*q(:))
770 : enddo
771 :
772 0 : r=sqrt(xx(1)**2+xx(2)**2+xx(3)**2)
773 :
774 0 : lam=sphere%lon; th=sphere%lat;
775 0 : sinlam=sin(lam); sinth=sin(th);
776 0 : coslam=cos(lam); costh=cos(th);
777 :
778 0 : D1(1,1)=-sinlam; D1(1,2)=coslam; D1(1,3)=0.0_r8;
779 0 : D1(2,1)=0.0_r8; D1(2,2)=0.0_r8; D1(2,3)=1.0_r8;
780 :
781 0 : D2(1,1)=(sinlam**2)*(costh**2)+sinth**2; D2(1,2)=-sinlam*coslam*(costh**2); D2(1,3)=-coslam*sinth*costh;
782 0 : D2(2,1)=-sinlam*coslam*(costh**2); D2(2,2)=(coslam**2)*(costh**2)+sinth**2; D2(2,3)=-sinlam*sinth*costh;
783 0 : D2(3,1)=-coslam*sinth; D2(3,2)=-sinlam*sinth; D2(3,3)=costh;
784 :
785 0 : dd(1,1)=-1+b; dd(1,2)=-1+a;
786 0 : dd(2,1)=1-b; dd(2,2)=-1-a;
787 0 : dd(3,1)=1+b; dd(3,2)=1+a;
788 0 : dd(4,1)=-1-b; dd(4,2)=1-a;
789 :
790 0 : dd=dd/4.0_r8
791 :
792 0 : do i=1,3
793 0 : do j=1,2
794 0 : D3(i,j)=sum(c(i,:)*dd(:,j))
795 : enddo
796 : enddo
797 :
798 0 : do i=1,3
799 0 : do j=1,2
800 0 : D4(i,j)=sum(D2(i,:)*D3(:,j))
801 : enddo
802 : enddo
803 :
804 0 : do i=1,2
805 0 : do j=1,2
806 0 : D(i,j)=sum(D1(i,:)*D4(:,j))
807 : enddo
808 : enddo
809 :
810 0 : D=D/r
811 0 : end subroutine dmap_elementlocal
812 :
813 :
814 :
815 :
816 :
817 : ! ========================================
818 : ! coreolis_init_atomic:
819 : !
820 : ! Initialize coreolis term ...
821 : !
822 : ! ========================================
823 :
824 21600 : subroutine coreolis_init_atomic(elem)
825 : use element_mod, only: element_t
826 : use dimensions_mod, only: np
827 : use physconst, only: omega
828 :
829 : type (element_t) :: elem
830 :
831 : ! Local variables
832 :
833 : integer :: i,j
834 : real (kind=r8) :: lat,lon,rangle
835 :
836 21600 : rangle = rotate_grid * PI / 180._r8
837 108000 : do j=1,np
838 453600 : do i=1,np
839 432000 : if ( rotate_grid /= 0) then
840 0 : lat = elem%spherep(i,j)%lat
841 0 : lon = elem%spherep(i,j)%lon
842 : elem%fcor(i,j)= 2*omega* &
843 0 : (-cos(lon)*cos(lat)*sin(rangle) + sin(lat)*cos(rangle))
844 : else
845 345600 : elem%fcor(i,j) = 2.0_r8*omega*SIN(elem%spherep(i,j)%lat)
846 : endif
847 : end do
848 : end do
849 :
850 21600 : end subroutine coreolis_init_atomic
851 :
852 : ! =========================================
853 : ! rotation_init_atomic:
854 : !
855 : ! Initialize cube rotation terms resulting
856 : ! from changing cube face coordinate systems
857 : !
858 : ! =========================================
859 :
860 :
861 10800 : subroutine rotation_init_atomic(elem, rot_type)
862 : use element_mod, only : element_t
863 : use dimensions_mod, only : np
864 : use control_mod, only : north, south, east, west, neast, seast, swest, nwest
865 :
866 : type (element_t) :: elem
867 : character(len=*) rot_type
868 :
869 : ! =======================================
870 : ! Local variables
871 : ! =======================================
872 :
873 : integer :: myface_no ! current element face number
874 : integer :: nbrface_no ! neighbor element face number
875 : integer :: inbr
876 : integer :: nrot,irot
877 : integer :: ii,i,j,k
878 : integer :: ir,jr
879 : integer :: start, cnt
880 :
881 : real (kind=r8) :: Dloc(2,2,np)
882 : real (kind=r8) :: Drem(2,2,np)
883 : real (kind=r8) :: x1,x2
884 :
885 :
886 10800 : myface_no = elem%vertex%face_number
887 :
888 10800 : nrot = 0
889 :
890 97200 : do inbr=1,8
891 86400 : cnt = elem%vertex%nbrs_ptr(inbr+1) - elem%vertex%nbrs_ptr(inbr)
892 : start = elem%vertex%nbrs_ptr(inbr)
893 :
894 183552 : do k = 0, cnt-1
895 86352 : nbrface_no = elem%vertex%nbrs_face(start+k)
896 172752 : if (myface_no /= nbrface_no) nrot=nrot+1
897 : end do
898 :
899 : end do
900 :
901 10800 : if(associated(elem%desc%rot)) then
902 0 : if (size(elem%desc%rot) > 0) then
903 : ! deallocate(elem%desc%rot)
904 0 : NULLIFY(elem%desc%rot)
905 : endif
906 : end if
907 :
908 : ! =====================================================
909 : ! If there are neighbors on other cube faces, allocate
910 : ! an array of rotation matrix structs.
911 : ! =====================================================
912 :
913 10800 : if (nrot > 0) then
914 8400 : allocate(elem%desc%rot(nrot))
915 1392 : elem%desc%use_rotation=1
916 1392 : irot=0
917 :
918 12528 : do inbr=1,8
919 11136 : cnt = elem%vertex%nbrs_ptr(inbr+1) - elem%vertex%nbrs_ptr(inbr)
920 : start = elem%vertex%nbrs_ptr(inbr)
921 :
922 23616 : do k= 0, cnt-1
923 :
924 11088 : nbrface_no = elem%vertex%nbrs_face(start+k)
925 : ! The cube edge (myface_no,nbrface_no) and inbr defines
926 : ! a unique rotation given by (D^-1) on myface_no x (D) on nbrface_no
927 :
928 22224 : if (myface_no /= nbrface_no .and. elem%vertex%nbrs(start+k) /= -1 ) then
929 4224 : irot=irot+1
930 :
931 4224 : if (inbr <= 4) then
932 1440 : allocate(elem%desc%rot(irot)%R(2,2,np)) ! edge
933 : else
934 2784 : allocate(elem%desc%rot(irot)%R(2,2,1 )) ! corner
935 : end if
936 : ! Initialize Dloc and Drem for no-rotation possibilities
937 21120 : Dloc(1,1,:) = 1.0_r8
938 21120 : Dloc(1,2,:) = 0.0_r8
939 21120 : Dloc(2,1,:) = 0.0_r8
940 21120 : Dloc(2,2,:) = 1.0_r8
941 21120 : Drem(1,1,:) = 1.0_r8
942 21120 : Drem(1,2,:) = 0.0_r8
943 21120 : Drem(2,1,:) = 0.0_r8
944 21120 : Drem(2,2,:) = 1.0_r8
945 :
946 : ! must compute Dloc on my face, Drem on neighbor face,
947 : ! for each point on edge or corner.
948 :
949 : ! ====================================
950 : ! Equatorial belt east/west neighbors
951 : ! ====================================
952 :
953 4224 : if (nbrface_no <= 4 .and. myface_no <= 4) then
954 :
955 : if (inbr == west) then
956 1200 : do j=1,np
957 960 : x1 = elem%cartp(1,j)%x
958 960 : x2 = elem%cartp(1,j)%y
959 : call Vmap(Dloc(1,1,j), x1,x2,myface_no)
960 1200 : call Vmap(Drem(1,1,j),-x1,x2,nbrface_no)
961 : end do
962 : else if (inbr == east) then
963 1200 : do j=1,np
964 960 : x1 = elem%cartp(np,j)%x
965 960 : x2 = elem%cartp(np,j)%y
966 960 : call Vmap(Dloc(1,1,j), x1,x2,myface_no)
967 1200 : call Vmap(Drem(1,1,j),-x1,x2,nbrface_no)
968 : end do
969 : else if (inbr == swest ) then
970 232 : x1 = elem%cartp(1,1)%x
971 232 : x2 = elem%cartp(1,1)%y
972 232 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
973 232 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
974 : else if (inbr == nwest ) then
975 232 : x1 = elem%cartp(1,np)%x
976 232 : x2 = elem%cartp(1,np)%y
977 232 : call Vmap(Dloc(1,1,1), x1,x2,myface_no)
978 232 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
979 : else if (inbr == seast ) then
980 232 : x1 = elem%cartp(np,1)%x
981 232 : x2 = elem%cartp(np,1)%y
982 232 : call Vmap(Dloc(1,1,1), x1,x2,myface_no)
983 232 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
984 : else if (inbr == neast ) then
985 232 : x1 = elem%cartp(np,np)%x
986 232 : x2 = elem%cartp(np,np)%y
987 232 : call Vmap(Dloc(1,1,1), x1,x2,myface_no)
988 232 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
989 : end if
990 :
991 : end if
992 :
993 : ! Northern Neighbors of Equatorial Belt
994 :
995 4224 : if ( myface_no <= 4 .and. nbrface_no == 6 ) then
996 704 : if (inbr == north) then
997 1200 : do i=1,np
998 960 : ir=np+1-i
999 960 : x1 = elem%cartp(i,np)%x
1000 960 : x2 = elem%cartp(i,np)%y
1001 960 : if ( myface_no == 1) then
1002 240 : call Vmap(Dloc(1,1,i), x1,x2,myface_no)
1003 240 : call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
1004 : end if
1005 1200 : if ( myface_no == 2) then
1006 240 : call Vmap(Dloc(1,1,i),x1,x2,myface_no)
1007 : call Vmap(Drem(1,1,i),x2,x1,nbrface_no)
1008 :
1009 : end if
1010 960 : if ( myface_no == 3) then
1011 240 : call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
1012 240 : call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
1013 : end if
1014 1200 : if ( myface_no == 4) then
1015 240 : call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
1016 240 : call Vmap(Drem(1,1,ir),-x2,-x1,nbrface_no)
1017 : end if
1018 : end do
1019 464 : else if (inbr == nwest) then
1020 232 : x1 = elem%cartp(1,np)%x
1021 232 : x2 = elem%cartp(1,np)%y
1022 232 : call Vmap(Dloc(1,1,1), x1,x2,myface_no)
1023 232 : if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1024 232 : if ( myface_no == 2) call Vmap(Drem(1,1,1),x2, x1,nbrface_no)
1025 232 : if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1026 232 : if ( myface_no == 4) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1027 232 : else if (inbr == neast) then
1028 232 : x1 = elem%cartp(np,np)%x
1029 232 : x2 = elem%cartp(np,np)%y
1030 232 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1031 232 : if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1032 232 : if ( myface_no == 2) call Vmap(Drem(1,1,1),x2, x1,nbrface_no)
1033 232 : if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1034 232 : if ( myface_no == 4) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1035 : end if
1036 :
1037 : end if
1038 :
1039 : ! Southern Neighbors of Equatorial Belt
1040 :
1041 4224 : if ( myface_no <= 4 .and. nbrface_no == 5 ) then
1042 704 : if (inbr == south) then
1043 1200 : do i=1,np
1044 960 : ir=np+1-i
1045 960 : x1 = elem%cartp(i,1)%x
1046 960 : x2 = elem%cartp(i,1)%y
1047 960 : if ( myface_no == 1) then
1048 240 : call Vmap(Dloc(1,1,i), x1, x2,myface_no)
1049 240 : call Vmap(Drem(1,1,i), x1,-x2,nbrface_no)
1050 : end if
1051 960 : if ( myface_no == 2) then
1052 240 : call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
1053 240 : call Vmap(Drem(1,1,ir),-x2,-x1,nbrface_no)
1054 : end if
1055 960 : if ( myface_no == 3) then
1056 240 : call Vmap(Dloc(1,1,ir), x1,x2,myface_no)
1057 240 : call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
1058 : end if
1059 1440 : if ( myface_no == 4) then
1060 240 : call Vmap(Dloc(1,1,i), x1,x2,myface_no)
1061 : call Vmap(Drem(1,1,i), x2,x1,nbrface_no)
1062 : end if
1063 : end do
1064 464 : else if (inbr == swest) then
1065 232 : x1 = elem%cartp(1,1)%x
1066 232 : x2 = elem%cartp(1,1)%y
1067 232 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1068 :
1069 :
1070 232 : if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1071 232 : if ( myface_no == 2) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1072 232 : if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1073 232 : if ( myface_no == 4) call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1074 :
1075 232 : else if (inbr == seast) then
1076 232 : x1 = elem%cartp(np,1)%x
1077 232 : x2 = elem%cartp(np,1)%y
1078 232 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1079 232 : if ( myface_no == 1) call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1080 232 : if ( myface_no == 2) call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1081 232 : if ( myface_no == 3) call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1082 232 : if ( myface_no == 4) call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1083 : end if
1084 :
1085 : end if
1086 :
1087 : ! Neighbors of Northern Capping Face Number 6
1088 :
1089 4224 : if ( myface_no == 6 ) then
1090 704 : if (nbrface_no == 1) then
1091 176 : if (inbr == south) then
1092 300 : do i=1,np
1093 240 : x1 = elem%cartp(i,1)%x
1094 240 : x2 = elem%cartp(i,1)%y
1095 240 : call Vmap(Dloc(1,1,i),x1,x2,myface_no)
1096 300 : call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
1097 : end do
1098 116 : else if (inbr == swest) then
1099 58 : x1 = elem%cartp(1,1)%x
1100 58 : x2 = elem%cartp(1,1)%y
1101 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1102 58 : call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1103 58 : else if (inbr == seast) then
1104 58 : x1 = elem%cartp(np,1)%x
1105 58 : x2 = elem%cartp(np,1)%y
1106 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1107 58 : call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1108 : end if
1109 528 : else if (nbrface_no == 2) then
1110 176 : if (inbr == east) then
1111 300 : do j=1,np
1112 240 : x1 = elem%cartp(np,j)%x
1113 240 : x2 = elem%cartp(np,j)%y
1114 240 : call Vmap(Dloc(1,1,j),x1,x2,myface_no)
1115 300 : call Vmap(Drem(1,1,j),x2,x1,nbrface_no)
1116 : end do
1117 116 : else if (inbr == seast) then
1118 58 : x1 = elem%cartp(np,1)%x
1119 58 : x2 = elem%cartp(np,1)%y
1120 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1121 58 : call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1122 58 : else if (inbr == neast) then
1123 58 : x1 = elem%cartp(np,np)%x
1124 58 : x2 = elem%cartp(np,np)%y
1125 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1126 58 : call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1127 : end if
1128 352 : else if (nbrface_no == 3) then
1129 176 : if (inbr == north) then
1130 300 : do i=1,np
1131 240 : ir =np+1-i
1132 240 : x1 = elem%cartp(i,np)%x
1133 240 : x2 = elem%cartp(i,np)%y
1134 240 : call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
1135 300 : call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
1136 : end do
1137 116 : else if (inbr == nwest) then
1138 58 : x1 = elem%cartp(1,np)%x
1139 58 : x2 = elem%cartp(1,np)%y
1140 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1141 58 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1142 58 : else if (inbr == neast) then
1143 58 : x1 = elem%cartp(np,np)%x
1144 58 : x2 = elem%cartp(np,np)%y
1145 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1146 58 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1147 : end if
1148 176 : else if (nbrface_no == 4) then
1149 176 : if (inbr == west) then
1150 300 : do j=1,np
1151 240 : jr=np+1-j
1152 240 : x1 = elem%cartp(1,j)%x
1153 240 : x2 = elem%cartp(1,j)%y
1154 240 : call Vmap(Dloc(1,1,jr), x1, x2,myface_no )
1155 300 : call Vmap(Drem(1,1,jr),-x2,-x1,nbrface_no)
1156 : end do
1157 116 : else if (inbr == swest) then
1158 58 : x1 = elem%cartp(1,1)%x
1159 58 : x2 = elem%cartp(1,1)%y
1160 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1161 58 : call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1162 58 : else if (inbr == nwest) then
1163 58 : x1 = elem%cartp(1,np)%x
1164 58 : x2 = elem%cartp(1,np)%y
1165 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1166 58 : call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1167 : end if
1168 : end if
1169 : end if
1170 :
1171 : ! Neighbors of South Capping Face Number 5
1172 :
1173 4224 : if ( myface_no == 5 ) then
1174 704 : if (nbrface_no == 1) then
1175 176 : if (inbr == north) then
1176 300 : do i=1,np
1177 240 : x1 = elem%cartp(i,np)%x
1178 240 : x2 = elem%cartp(i,np)%y
1179 240 : call Vmap(Dloc(1,1,i),x1,x2,myface_no)
1180 300 : call Vmap(Drem(1,1,i),x1,-x2,nbrface_no)
1181 : end do
1182 116 : else if (inbr == nwest) then
1183 58 : x1 = elem%cartp(1,np)%x
1184 58 : x2 = elem%cartp(1,np)%y
1185 58 : call Vmap(Dloc(:,:,1),x1,x2,myface_no)
1186 58 : call Vmap(Drem(:,:,1),x1,-x2,nbrface_no)
1187 58 : else if (inbr == neast) then
1188 58 : x1 = elem%cartp(np,np)%x
1189 58 : x2 = elem%cartp(np,np)%y
1190 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1191 58 : call Vmap(Drem(1,1,1),x1,-x2,nbrface_no)
1192 : end if
1193 528 : else if (nbrface_no == 2) then
1194 176 : if (inbr == east) then
1195 300 : do j=1,np
1196 240 : jr=np+1-j
1197 240 : x1 = elem%cartp(np,j)%x
1198 240 : x2 = elem%cartp(np,j)%y
1199 240 : call Vmap(Dloc(1,1,jr),x1, x2,myface_no)
1200 300 : call Vmap(Drem(1,1,jr),-x2,-x1,nbrface_no)
1201 : end do
1202 116 : else if (inbr == seast) then
1203 58 : x1 = elem%cartp(np,1)%x
1204 58 : x2 = elem%cartp(np,1)%y
1205 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1206 58 : call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1207 58 : else if (inbr == neast) then
1208 58 : x1 = elem%cartp(np,np)%x
1209 58 : x2 = elem%cartp(np,np)%y
1210 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1211 58 : call Vmap(Drem(1,1,1),-x2,-x1,nbrface_no)
1212 : end if
1213 352 : else if (nbrface_no == 3) then
1214 176 : if (inbr == south) then
1215 300 : do i=1,np
1216 240 : ir=np+1-i
1217 240 : x1 = elem%cartp(i,1)%x
1218 240 : x2 = elem%cartp(i,1)%y
1219 240 : call Vmap(Dloc(1,1,ir),x1,x2,myface_no)
1220 300 : call Vmap(Drem(1,1,ir),-x1,x2,nbrface_no)
1221 : end do
1222 116 : else if (inbr == swest) then
1223 58 : x1 = elem%cartp(1,1)%x
1224 58 : x2 = elem%cartp(1,1)%y
1225 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1226 58 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1227 58 : else if (inbr == seast) then
1228 58 : x1 = elem%cartp(np,1)%x
1229 58 : x2 = elem%cartp(np,1)%y
1230 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1231 58 : call Vmap(Drem(1,1,1),-x1,x2,nbrface_no)
1232 : end if
1233 176 : else if (nbrface_no == 4) then
1234 176 : if (inbr == west) then
1235 300 : do j=1,np
1236 240 : x1 = elem%cartp(1,j)%x
1237 240 : x2 = elem%cartp(1,j)%y
1238 : call Vmap(Dloc(1,1,j),x1,x2,myface_no)
1239 300 : call Vmap(Drem(1,1,j),x2,x1,nbrface_no)
1240 : end do
1241 116 : else if (inbr == swest) then
1242 58 : x1 = elem%cartp(1,1)%x
1243 58 : x2 = elem%cartp(1,1)%y
1244 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1245 58 : call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1246 58 : else if (inbr == nwest) then
1247 58 : x1 = elem%cartp(1,np)%x
1248 58 : x2 = elem%cartp(1,np)%y
1249 58 : call Vmap(Dloc(1,1,1),x1,x2,myface_no)
1250 58 : call Vmap(Drem(1,1,1),x2,x1,nbrface_no)
1251 : end if
1252 : end if
1253 : end if
1254 :
1255 4224 : elem%desc%rot(irot)%nbr = inbr
1256 4224 : if (rot_type == "covariant") then
1257 0 : do i=1,SIZE(elem%desc%rot(irot)%R(:,:,:),3)
1258 0 : elem%desc%rot(irot)%R(:,:,i)=covariant_rot(Dloc(:,:,i),Drem(:,:,i))
1259 : end do
1260 4224 : else if (rot_type == "contravariant") then
1261 12768 : do i=1,SIZE(elem%desc%rot(irot)%R(:,:,:),3)
1262 64032 : elem%desc%rot(irot)%R(:,:,i)=contravariant_rot(Dloc(:,:,i),Drem(:,:,i))
1263 : end do
1264 : end if
1265 :
1266 : end if ! end of a unique rotation
1267 : end do !k loop over neighbors in that direction
1268 : end do !inbr loop
1269 : end if !nrot > 0
1270 :
1271 10800 : end subroutine rotation_init_atomic
1272 :
1273 :
1274 21600 : subroutine set_corner_coordinates(elem)
1275 : use element_mod, only : element_t
1276 : use dimensions_mod, only : ne
1277 :
1278 : type (element_t) :: elem
1279 :
1280 : ! Local variables
1281 : integer i,ie,je,face_no,nn
1282 : real (kind=r8) :: dx,dy, startx, starty
1283 :
1284 10800 : if (0==ne) call endrun('Error in set_corner_coordinates: ne is zero')
1285 :
1286 : ! ========================================
1287 : ! compute cube face coordinates of element
1288 : ! =========================================
1289 :
1290 10800 : call convert_gbl_index(elem%vertex%number,ie,je,face_no)
1291 :
1292 10800 : elem%vertex%face_number = face_no
1293 10800 : dx = (cube_xend-cube_xstart)/ne
1294 10800 : dy = (cube_yend-cube_ystart)/ne
1295 :
1296 10800 : startx = cube_xstart+ie*dx
1297 10800 : starty = cube_ystart+je*dy
1298 :
1299 10800 : elem%corners(1)%x = startx
1300 10800 : elem%corners(1)%y = starty
1301 10800 : elem%corners(2)%x = startx+dx
1302 10800 : elem%corners(2)%y = starty
1303 10800 : elem%corners(3)%x = startx+dx
1304 10800 : elem%corners(3)%y = starty+dy
1305 10800 : elem%corners(4)%x = startx
1306 10800 : elem%corners(4)%y = starty+dy
1307 :
1308 : #if 0
1309 : do i=1,4
1310 : elem%node_multiplicity(i) = 4
1311 : end do
1312 : ie = ie + 1
1313 : je = je + 1
1314 : if (ie == 1 .and. je == 1) then
1315 : elem%node_multiplicity(1) = 3
1316 : else if (ie == ne .and. je == 1) then
1317 : elem%node_multiplicity(2) = 3
1318 : else if (ie == ne .and. je == ne) then
1319 : elem%node_multiplicity(3) = 3
1320 : else if (ie == 1 .and. je == ne) then
1321 : elem%node_multiplicity(4) = 3
1322 : end if
1323 : #endif
1324 10800 : end subroutine set_corner_coordinates
1325 :
1326 :
1327 1536 : subroutine assign_node_numbers_to_elem(elements, GridVertex)
1328 : use dimensions_mod, only : ne
1329 : use element_mod, only : element_t
1330 : use control_mod, only : north, south, east, west, neast, seast, swest, nwest
1331 : use gridgraph_mod, only : GridVertex_t
1332 : implicit none
1333 : type (element_t), intent(inout) :: elements(:)
1334 : type (GridVertex_t), intent(in) :: GridVertex(:)
1335 :
1336 : type (GridVertex_t) :: vertex
1337 3072 : integer :: connectivity(6*ne*ne, 4)
1338 : integer :: nn(4), en(4)
1339 : integer el, i, n, direction
1340 : integer current_node_num, tot_ne
1341 : integer :: start, cnt
1342 :
1343 1536 : current_node_num = 0
1344 1536 : tot_ne = 6*ne*ne
1345 :
1346 0 : if (0==ne) call endrun('Error in assign_node_numbers_to_elem: ne is zero')
1347 1536 : if (tot_ne /= SIZE(GridVertex)) call endrun('Error in assign_node_numbers_to_elem: GridVertex not correct length')
1348 :
1349 33185280 : connectivity = 0
1350 :
1351 8295936 : do el = 1,tot_ne
1352 8294400 : vertex = GridVertex(el)
1353 8294400 : en = 0
1354 74649600 : do direction = 1,8
1355 66355200 : cnt = vertex%nbrs_ptr(direction+1) - vertex%nbrs_ptr(direction)
1356 : start = vertex%nbrs_ptr(direction)
1357 :
1358 140967936 : do i=0, cnt-1
1359 66318336 : n = vertex%nbrs(start+i)
1360 132673536 : if (n /= -1) then
1361 331591680 : nn = connectivity(n,:)
1362 8294400 : select case (direction)
1363 : case (north)
1364 8294400 : if (nn(1)/=0) en(4) = nn(1)
1365 8294400 : if (nn(2)/=0) en(3) = nn(2)
1366 : case (south)
1367 8294400 : if (nn(4)/=0) en(1) = nn(4)
1368 8294400 : if (nn(3)/=0) en(2) = nn(3)
1369 : case (east)
1370 8294400 : if (nn(1)/=0) en(2) = nn(1)
1371 8294400 : if (nn(4)/=0) en(3) = nn(4)
1372 : case (west)
1373 8294400 : if (nn(2)/=0) en(1) = nn(2)
1374 8294400 : if (nn(3)/=0) en(4) = nn(3)
1375 : case (neast)
1376 8285184 : if (nn(1)/=0) en(3) = nn(1)
1377 : case (seast)
1378 8285184 : if (nn(4)/=0) en(2) = nn(4)
1379 : case (swest)
1380 8285184 : if (nn(3)/=0) en(1) = nn(3)
1381 : case (nwest)
1382 66318336 : if (nn(2)/=0) en(4) = nn(2)
1383 : end select
1384 : end if
1385 : end do
1386 : end do !direction
1387 :
1388 41472000 : do i=1,4
1389 41472000 : if (en(i) == 0) then
1390 8297472 : current_node_num = current_node_num + 1
1391 8297472 : en(i) = current_node_num
1392 : end if
1393 : end do
1394 41473536 : connectivity(el,:) = en
1395 : end do
1396 :
1397 1536 : if (current_node_num /= (6*ne*ne+2)) then
1398 0 : call endrun('Error in assignment of node numbers: Failed Euler test')
1399 : end if
1400 : ! do el = 1,SIZE(elements)
1401 : ! elements(el)%node_numbers = connectivity(elements(el)%vertex%number, :)
1402 : ! end do
1403 1536 : end subroutine assign_node_numbers_to_elem
1404 :
1405 :
1406 : ! ================================================
1407 : ! convert_gbl_index:
1408 : !
1409 : ! Convert global element index to cube index
1410 : ! ================================================
1411 :
1412 10800 : subroutine convert_gbl_index(number,ie,je,face_no)
1413 : use dimensions_mod, only : ne
1414 : integer, intent(in) :: number
1415 : integer, intent(out) :: ie,je,face_no
1416 :
1417 10800 : if (0==ne) call endrun('Error in cube_mod:convert_gbl_index: ne is zero')
1418 :
1419 : ! inverse of the function: number = 1 + ie + ne*je + ne*ne*(face_no-1)
1420 10800 : face_no=((number-1)/(ne*ne))+1
1421 10800 : ie=MODULO(number-1,ne)
1422 10800 : je=(number-1)/ne - (face_no-1)*ne
1423 :
1424 10800 : end subroutine convert_gbl_index
1425 :
1426 1536 : subroutine CubeTopology(GridEdge, GridVertex)
1427 : use gridgraph_mod, only: GridEdge_t, GridVertex_t, initgridedge
1428 : use gridgraph_mod, only: allocate_gridvertex_nbrs, deallocate_gridvertex_nbrs
1429 : use dimensions_mod, only: np, ne
1430 : use spacecurve_mod, only: IsFactorable, genspacecurve
1431 : use control_mod, only: north, south, east, west, neast, seast, swest, nwest
1432 : !-----------------------
1433 :
1434 : ! Since GridVertex fields must be allocated before calling this, it
1435 : ! must be intent(inout).
1436 : !og: is 'target' here necessary?
1437 : !GridEdge : changed its 'out' attribute to 'inout'
1438 : type (GridEdge_t), intent(inout),target :: GridEdge(:)
1439 : type (GridVertex_t), intent(inout),target :: GridVertex(:)
1440 :
1441 :
1442 1536 : integer,allocatable :: Mesh(:,:)
1443 1536 : integer,allocatable :: Mesh2(:,:),Mesh2_map(:,:,:),sfcij(:,:)
1444 1536 : type (GridVertex_t),allocatable :: GridElem(:,:,:)
1445 : integer :: i,j,k,ll,number,irev,ne2,i2,j2,sfc_index
1446 : integer :: EdgeWgtP,CornerWgt
1447 : integer :: ielem, nedge
1448 : integer :: offset, ierr, loc
1449 1536 : logical, allocatable :: nbrs_used(:,:,:,:)
1450 :
1451 :
1452 1536 : if (0==ne) call endrun('Error in CubeTopology: ne is zero')
1453 :
1454 8587776 : allocate(GridElem(ne,ne,nfaces),stat=ierr)
1455 10752 : do k = 1, nfaces
1456 287232 : do j = 1, ne
1457 8580096 : do i = 1, ne
1458 8570880 : call allocate_gridvertex_nbrs(GridElem(i,j,k))
1459 : end do
1460 : end do
1461 : end do
1462 :
1463 1536 : if(ierr/=0) then
1464 0 : call endrun('error in allocation of GridElem structure')
1465 : end if
1466 :
1467 9216 : allocate(nbrs_used(ne,ne,nfaces,8))
1468 68654592 : nbrs_used = .false.
1469 :
1470 :
1471 : number=1
1472 : EdgeWgtP = np
1473 : CornerWgt = 1
1474 10752 : do k=1,nfaces
1475 287232 : do j=1,ne
1476 8580096 : do i=1,ne
1477 : ! ====================================
1478 : ! Number elements
1479 : ! ====================================
1480 74649600 : GridElem(i,j,k)%nbrs(:)=0
1481 74649600 : GridElem(i,j,k)%nbrs_wgt(:)=0
1482 82944000 : GridElem(i,j,k)%nbrs_ptr(:)=0
1483 74649600 : GridElem(i,j,k)%nbrs_wgt_ghost(:)=1 ! always this value
1484 8294400 : GridElem(i,j,k)%SpaceCurve=0
1485 8294400 : GridElem(i,j,k)%number=number
1486 8570880 : number=number+1
1487 :
1488 : end do
1489 : end do
1490 : end do
1491 :
1492 6144 : allocate(Mesh(ne,ne))
1493 1536 : if(IsFactorable(ne)) then
1494 1536 : call GenspaceCurve(Mesh)
1495 : else
1496 : ! find the smallest ne2 which is a power of 2 and ne2>ne
1497 0 : ne2 = 2**ceiling(log(real(ne)) / log(2.0_r8))
1498 0 : if (ne2 < ne) then
1499 0 : call endrun('Fatal SFC error')
1500 : end if
1501 :
1502 0 : allocate(Mesh2(ne2,ne2))
1503 0 : allocate(Mesh2_map(ne2,ne2,2))
1504 0 : allocate(sfcij(0:ne2*ne2,2))
1505 :
1506 0 : call GenspaceCurve(Mesh2) ! SFC partition for ne2
1507 :
1508 : ! associate every element on the ne x ne mesh (Mesh)
1509 : ! with its closest element on the ne2 x ne2 mesh (Mesh2)
1510 : ! Store this as a map from Mesh2 -> Mesh in Mesh2_map.
1511 : ! elements in Mesh2 which are not mapped get assigned a value of 0
1512 0 : Mesh2_map=0
1513 0 : do j=1,ne
1514 0 : do i=1,ne
1515 : ! map this element to an (i2,j2) element
1516 : ! [ (i-.5)/ne , (j-.5)/ne ] = [ (i2-.5)/ne2 , (j2-.5)/ne2 ]
1517 0 : i2=nint( ((i-0.5_r8)/ne)*ne2 + 0.5_r8 )
1518 0 : j2=nint( ((j-0.5_r8)/ne)*ne2 + 0.5_r8 )
1519 0 : if (i2<1) i2=1
1520 0 : if (i2>ne2) i2=ne2
1521 0 : if (j2<1) j2=1
1522 0 : if (j2>ne2) j2=ne2
1523 0 : Mesh2_map(i2,j2,1)=i
1524 0 : Mesh2_map(i2,j2,2)=j
1525 : enddo
1526 : enddo
1527 :
1528 : ! create a reverse index array for Mesh2
1529 : ! k = Mesh2(i,j)
1530 : ! (i,j) = (sfcij(k,1),sfci(k,2))
1531 0 : do j=1,ne2
1532 0 : do i=1,ne2
1533 0 : k=Mesh2(i,j)
1534 0 : sfcij(k,1)=i
1535 0 : sfcij(k,2)=j
1536 : enddo
1537 : enddo
1538 :
1539 : ! generate a SFC for Mesh with the same ordering as the
1540 : ! elements in Mesh2 which map to Mesh.
1541 : sfc_index=0
1542 0 : do k=0,ne2*ne2-1
1543 0 : i2=sfcij(k,1)
1544 0 : j2=sfcij(k,2)
1545 0 : i=Mesh2_map(i2,j2,1)
1546 0 : j=Mesh2_map(i2,j2,2)
1547 0 : if (i/=0) then
1548 : ! (i2,j2) element maps to (i,j) element
1549 0 : Mesh(i,j)=sfc_index
1550 0 : sfc_index=sfc_index+1
1551 : endif
1552 : enddo
1553 :
1554 0 : deallocate(Mesh2)
1555 0 : deallocate(Mesh2_map)
1556 0 : deallocate(sfcij)
1557 : endif
1558 :
1559 : ! -------------------------------------------
1560 : ! Setup the space-filling curve for face 1
1561 : ! -------------------------------------------
1562 1536 : offset=0
1563 47616 : do j=1,ne
1564 1430016 : do i=1,ne
1565 1428480 : GridElem(i,j,1)%SpaceCurve = offset + Mesh(i,ne-j+1)
1566 : enddo
1567 : enddo
1568 :
1569 : ! -------------------------------------------
1570 : ! Setup the space-filling curve for face 2
1571 : ! -------------------------------------------
1572 1536 : offset = offset + ne*ne
1573 47616 : do j=1,ne
1574 1430016 : do i=1,ne
1575 1428480 : GridElem(i,j,2)%SpaceCurve = offset + Mesh(i,ne-j+1)
1576 : enddo
1577 : enddo
1578 :
1579 : ! -------------------------------------------
1580 : ! Setup the space-filling curve for face 6
1581 : ! -------------------------------------------
1582 1536 : offset = offset + ne*ne
1583 47616 : do j=1,ne
1584 1430016 : do i=1,ne
1585 1428480 : GridElem(i,j,6)%SpaceCurve = offset + Mesh(ne-i+1,ne-j+1)
1586 : enddo
1587 : enddo
1588 :
1589 : ! -------------------------------------------
1590 : ! Setup the space-filling curve for face 4
1591 : ! -------------------------------------------
1592 1536 : offset = offset + ne*ne
1593 47616 : do j=1,ne
1594 1430016 : do i=1,ne
1595 1428480 : GridElem(i,j,4)%SpaceCurve = offset + Mesh(ne-j+1,i)
1596 : enddo
1597 : enddo
1598 :
1599 : ! -------------------------------------------
1600 : ! Setup the space-filling curve for face 5
1601 : ! -------------------------------------------
1602 1536 : offset = offset + ne*ne
1603 47616 : do j=1,ne
1604 1430016 : do i=1,ne
1605 1428480 : GridElem(i,j,5)%SpaceCurve = offset + Mesh(i,j)
1606 : enddo
1607 : enddo
1608 :
1609 :
1610 : ! -------------------------------------------
1611 : ! Setup the space-filling curve for face 3
1612 : ! -------------------------------------------
1613 1536 : offset = offset + ne*ne
1614 47616 : do j=1,ne
1615 1430016 : do i=1,ne
1616 1428480 : GridElem(i,j,3)%SpaceCurve = offset + Mesh(i,j)
1617 : enddo
1618 : enddo
1619 :
1620 : ! ==================
1621 : ! face interiors
1622 : ! ==================
1623 10752 : do k=1,6
1624 : ! setup SOUTH, WEST, SW neighbors
1625 276480 : do j=2,ne
1626 8027136 : do i=2,ne
1627 7750656 : nbrs_used(i,j,k,west) = .true.
1628 7750656 : nbrs_used(i,j,k,south) = .true.
1629 7750656 : nbrs_used(i,j,k,swest) = .true.
1630 :
1631 :
1632 7750656 : GridElem(i,j,k)%nbrs(west) = GridElem(i-1,j,k)%number
1633 7750656 : GridElem(i,j,k)%nbrs_face(west) = k
1634 7750656 : GridElem(i,j,k)%nbrs_wgt(west) = EdgeWgtP
1635 7750656 : GridElem(i,j,k)%nbrs(south) = GridElem(i,j-1,k)%number
1636 7750656 : GridElem(i,j,k)%nbrs_face(south) = k
1637 7750656 : GridElem(i,j,k)%nbrs_wgt(south) = EdgeWgtP
1638 7750656 : GridElem(i,j,k)%nbrs(swest) = GridElem(i-1,j-1,k)%number
1639 7750656 : GridElem(i,j,k)%nbrs_face(swest) = k
1640 8017920 : GridElem(i,j,k)%nbrs_wgt(swest) = CornerWgt
1641 : end do
1642 : end do
1643 :
1644 : ! setup EAST, NORTH, NE neighbors
1645 276480 : do j=1,ne-1
1646 8027136 : do i=1,ne-1
1647 7750656 : nbrs_used(i,j,k,east) = .true.
1648 7750656 : nbrs_used(i,j,k,north) = .true.
1649 7750656 : nbrs_used(i,j,k,neast) = .true.
1650 :
1651 7750656 : GridElem(i,j,k)%nbrs(east) = GridElem(i+1,j,k)%number
1652 7750656 : GridElem(i,j,k)%nbrs_face(east) = k
1653 7750656 : GridElem(i,j,k)%nbrs_wgt(east) = EdgeWgtP
1654 7750656 : GridElem(i,j,k)%nbrs(north) = GridElem(i,j+1,k)%number
1655 7750656 : GridElem(i,j,k)%nbrs_face(north) = k
1656 7750656 : GridElem(i,j,k)%nbrs_wgt(north) = EdgeWgtP
1657 7750656 : GridElem(i,j,k)%nbrs(neast) = GridElem(i+1,j+1,k)%number
1658 7750656 : GridElem(i,j,k)%nbrs_face(neast) = k
1659 8017920 : GridElem(i,j,k)%nbrs_wgt(neast) = CornerWgt
1660 : end do
1661 : end do
1662 :
1663 : ! Setup the remaining SOUTH, EAST, and SE neighbors
1664 276480 : do j=2,ne
1665 8027136 : do i=1,ne-1
1666 7750656 : nbrs_used(i,j,k,south) = .true.
1667 7750656 : nbrs_used(i,j,k,east) = .true.
1668 7750656 : nbrs_used(i,j,k,seast) = .true.
1669 :
1670 :
1671 :
1672 7750656 : GridElem(i,j,k)%nbrs(south) = GridElem(i,j-1,k)%number
1673 7750656 : GridElem(i,j,k)%nbrs_face(south) = k
1674 7750656 : GridElem(i,j,k)%nbrs_wgt(south) = EdgeWgtP
1675 7750656 : GridElem(i,j,k)%nbrs(east) = GridElem(i+1,j,k)%number
1676 7750656 : GridElem(i,j,k)%nbrs_face(east) = k
1677 7750656 : GridElem(i,j,k)%nbrs_wgt(east) = EdgeWgtP
1678 7750656 : GridElem(i,j,k)%nbrs(seast) = GridElem(i+1,j-1,k)%number
1679 7750656 : GridElem(i,j,k)%nbrs_face(seast) = k
1680 8017920 : GridElem(i,j,k)%nbrs_wgt(seast) = CornerWgt
1681 : enddo
1682 : enddo
1683 :
1684 : ! Setup the remaining NORTH, WEST, and NW neighbors
1685 278016 : do j=1,ne-1
1686 8027136 : do i=2,ne
1687 7750656 : nbrs_used(i,j,k,north) = .true.
1688 7750656 : nbrs_used(i,j,k,west) = .true.
1689 7750656 : nbrs_used(i,j,k,nwest) = .true.
1690 :
1691 :
1692 :
1693 7750656 : GridElem(i,j,k)%nbrs(north) = GridElem(i,j+1,k)%number
1694 7750656 : GridElem(i,j,k)%nbrs_face(north) = k
1695 7750656 : GridElem(i,j,k)%nbrs_wgt(north) = EdgeWgtP
1696 7750656 : GridElem(i,j,k)%nbrs(west) = GridElem(i-1,j,k)%number
1697 7750656 : GridElem(i,j,k)%nbrs_face(west) = k
1698 7750656 : GridElem(i,j,k)%nbrs_wgt(west) = EdgeWgtP
1699 7750656 : GridElem(i,j,k)%nbrs(nwest) = GridElem(i-1,j+1,k)%number
1700 7750656 : GridElem(i,j,k)%nbrs_face(nwest) = k
1701 8017920 : GridElem(i,j,k)%nbrs_wgt(nwest) = CornerWgt
1702 : enddo
1703 : enddo
1704 : end do
1705 :
1706 : ! ======================
1707 : ! west/east "belt" edges
1708 : ! ======================
1709 :
1710 7680 : do k=1,4
1711 192000 : do j=1,ne
1712 184320 : nbrs_used(1,j,k,west) = .true.
1713 184320 : nbrs_used(ne,j,k,east) = .true.
1714 :
1715 :
1716 184320 : GridElem(1 ,j,k)%nbrs(west) = GridElem(ne,j,MODULO(2+k,4)+1)%number
1717 184320 : GridElem(1 ,j,k)%nbrs_face(west) = MODULO(2+k,4)+1
1718 184320 : GridElem(1 ,j,k)%nbrs_wgt(west) = EdgeWgtP
1719 184320 : GridElem(ne,j,k)%nbrs(east) = GridElem(1 ,j,MODULO(k ,4)+1)%number
1720 184320 : GridElem(ne,j,k)%nbrs_face(east) = MODULO(k ,4)+1
1721 184320 : GridElem(ne,j,k)%nbrs_wgt(east) = EdgeWgtP
1722 :
1723 : ! Special rules for corner 'edges'
1724 184320 : if( j /= 1) then
1725 178176 : nbrs_used(1,j,k,swest) = .true.
1726 178176 : nbrs_used(ne,j,k,seast) = .true.
1727 :
1728 :
1729 178176 : GridElem(1 ,j,k)%nbrs(swest) = GridElem(ne,j-1,MODULO(2+k,4)+1)%number
1730 178176 : GridElem(1 ,j,k)%nbrs_face(swest) = MODULO(2+k,4)+1
1731 178176 : GridElem(1 ,j,k)%nbrs_wgt(swest) = CornerWgt
1732 178176 : GridElem(ne,j,k)%nbrs(seast) = GridElem(1 ,j-1,MODULO(k ,4)+1)%number
1733 178176 : GridElem(ne,j,k)%nbrs_face(seast) = MODULO(k ,4)+1
1734 178176 : GridElem(ne,j,k)%nbrs_wgt(seast) = CornerWgt
1735 : endif
1736 190464 : if( j /= ne) then
1737 178176 : nbrs_used(1,j,k,nwest) = .true.
1738 178176 : nbrs_used(ne,j,k,neast) = .true.
1739 :
1740 :
1741 178176 : GridElem(1 ,j,k)%nbrs(nwest) = GridElem(ne,j+1,MODULO(2+k,4)+1)%number
1742 178176 : GridElem(1 ,j,k)%nbrs_face(nwest) = MODULO(2+k,4)+1
1743 178176 : GridElem(1 ,j,k)%nbrs_wgt(nwest) = CornerWgt
1744 178176 : GridElem(ne,j,k)%nbrs(neast) = GridElem(1 ,j+1,MODULO(k ,4)+1)%number
1745 178176 : GridElem(ne,j,k)%nbrs_face(neast) = MODULO(k ,4)+1
1746 178176 : GridElem(ne,j,k)%nbrs_wgt(neast) = CornerWgt
1747 : endif
1748 : end do
1749 : end do
1750 :
1751 :
1752 : ! ==================================
1753 : ! south edge of 1 / north edge of 5
1754 : ! ==================================
1755 :
1756 47616 : do i=1,ne
1757 46080 : nbrs_used(i,1,1,south) = .true.
1758 46080 : nbrs_used(i,ne,5,north) = .true.
1759 :
1760 46080 : GridElem(i,1 ,1)%nbrs(south) = GridElem(i,ne,5)%number
1761 46080 : GridElem(i,1 ,1)%nbrs_face(south) = 5
1762 46080 : GridElem(i,1 ,1)%nbrs_wgt(south) = EdgeWgtP
1763 46080 : GridElem(i,ne,5)%nbrs(north) = GridElem(i,1 ,1)%number
1764 46080 : GridElem(i,ne,5)%nbrs_face(north) = 1
1765 46080 : GridElem(i,ne,5)%nbrs_wgt(north) = EdgeWgtP
1766 :
1767 : ! Special rules for corner 'edges'
1768 46080 : if( i /= 1) then
1769 44544 : nbrs_used(i,1,1,swest) = .true.
1770 44544 : nbrs_used(i,ne,5,nwest) = .true.
1771 :
1772 44544 : GridElem(i,1 ,1)%nbrs(swest) = GridElem(i-1,ne,5)%number
1773 44544 : GridElem(i,1 ,1)%nbrs_face(swest) = 5
1774 44544 : GridElem(i,1 ,1)%nbrs_wgt(swest) = CornerWgt
1775 44544 : GridElem(i,ne,5)%nbrs(nwest) = GridElem(i-1,1 ,1)%number
1776 44544 : GridElem(i,ne,5)%nbrs_face(nwest) = 1
1777 44544 : GridElem(i,ne,5)%nbrs_wgt(nwest) = CornerWgt
1778 : endif
1779 47616 : if( i /= ne) then
1780 44544 : nbrs_used(i,1,1,seast) = .true.
1781 44544 : nbrs_used(i,ne,5,neast) = .true.
1782 :
1783 44544 : GridElem(i,1 ,1)%nbrs(seast) = GridElem(i+1,ne,5)%number
1784 44544 : GridElem(i,1 ,1)%nbrs_face(seast) = 5
1785 44544 : GridElem(i,1 ,1)%nbrs_wgt(seast) = CornerWgt
1786 44544 : GridElem(i,ne,5)%nbrs(neast) = GridElem(i+1,1 ,1)%number
1787 44544 : GridElem(i,ne,5)%nbrs_face(neast) = 1
1788 44544 : GridElem(i,ne,5)%nbrs_wgt(neast) = CornerWgt
1789 : endif
1790 :
1791 : end do
1792 :
1793 : ! ==================================
1794 : ! south edge of 2 / east edge of 5
1795 : ! ==================================
1796 :
1797 47616 : do i=1,ne
1798 46080 : irev=ne+1-i
1799 46080 : nbrs_used(i,1,2,south) = .true.
1800 46080 : nbrs_used(ne,i,5,east) = .true.
1801 :
1802 :
1803 46080 : GridElem(i,1 ,2)%nbrs(south) = GridElem(ne,irev,5)%number
1804 46080 : GridElem(i,1 ,2)%nbrs_face(south) = 5
1805 46080 : GridElem(i,1 ,2)%nbrs_wgt(south) = EdgeWgtP
1806 46080 : GridElem(ne,i,5)%nbrs(east) = GridElem(irev,1 ,2)%number
1807 46080 : GridElem(ne,i,5)%nbrs_face(east) = 2
1808 46080 : GridElem(ne,i,5)%nbrs_wgt(east) = EdgeWgtP
1809 :
1810 : ! Special rules for corner 'edges'
1811 46080 : if( i /= 1) then
1812 44544 : nbrs_used(i,1,2,swest) = .true.
1813 44544 : nbrs_used(ne,i,5,seast) = .true.
1814 :
1815 :
1816 44544 : GridElem(i,1 ,2)%nbrs(swest) = GridElem(ne,irev+1,5)%number
1817 44544 : GridElem(i,1 ,2)%nbrs_face(swest) = 5
1818 44544 : GridElem(i,1 ,2)%nbrs_wgt(swest) = CornerWgt
1819 44544 : GridElem(ne,i,5)%nbrs(seast) = GridElem(irev+1,1 ,2)%number
1820 44544 : GridElem(ne,i,5)%nbrs_face(seast) = 2
1821 44544 : GridElem(ne,i,5)%nbrs_wgt(seast) = CornerWgt
1822 : endif
1823 47616 : if(i /= ne) then
1824 44544 : nbrs_used(i,1,2,seast) = .true.
1825 44544 : nbrs_used(ne,i,5,neast) = .true.
1826 :
1827 :
1828 44544 : GridElem(i,1 ,2)%nbrs(seast) = GridElem(ne,irev-1,5)%number
1829 44544 : GridElem(i,1 ,2)%nbrs_face(seast) = 5
1830 44544 : GridElem(i,1 ,2)%nbrs_wgt(seast) = CornerWgt
1831 44544 : GridElem(ne,i,5)%nbrs(neast) = GridElem(irev-1,1 ,2)%number
1832 44544 : GridElem(ne,i,5)%nbrs_face(neast) = 2
1833 44544 : GridElem(ne,i,5)%nbrs_wgt(neast) = CornerWgt
1834 : endif
1835 : enddo
1836 : ! ==================================
1837 : ! south edge of 3 / south edge of 5
1838 : ! ==================================
1839 :
1840 47616 : do i=1,ne
1841 46080 : irev=ne+1-i
1842 46080 : nbrs_used(i,1,3,south) = .true.
1843 46080 : nbrs_used(i,1,5,south) = .true.
1844 :
1845 46080 : GridElem(i,1,3)%nbrs(south) = GridElem(irev,1,5)%number
1846 46080 : GridElem(i,1,3)%nbrs_face(south) = 5
1847 46080 : GridElem(i,1,3)%nbrs_wgt(south) = EdgeWgtP
1848 46080 : GridElem(i,1,5)%nbrs(south) = GridElem(irev,1,3)%number
1849 46080 : GridElem(i,1,5)%nbrs_face(south) = 3
1850 46080 : GridElem(i,1,5)%nbrs_wgt(south) = EdgeWgtP
1851 :
1852 : ! Special rules for corner 'edges'
1853 46080 : if( i /= 1) then
1854 44544 : nbrs_used(i,1,3,swest) = .true.
1855 44544 : nbrs_used(i,1,5,swest) = .true.
1856 :
1857 :
1858 44544 : GridElem(i,1,3)%nbrs(swest) = GridElem(irev+1,1,5)%number
1859 44544 : GridElem(i,1,3)%nbrs_face(swest) = 5
1860 44544 : GridElem(i,1,3)%nbrs_wgt(swest) = CornerWgt
1861 44544 : GridElem(i,1,5)%nbrs(swest) = GridElem(irev+1,1,3)%number
1862 44544 : GridElem(i,1,5)%nbrs_face(swest) = 3
1863 44544 : GridElem(i,1,5)%nbrs_wgt(swest) = CornerWgt
1864 : endif
1865 47616 : if(i /= ne) then
1866 44544 : nbrs_used(i,1,3,seast) = .true.
1867 44544 : nbrs_used(i,1,5,seast) = .true.
1868 :
1869 44544 : GridElem(i,1,3)%nbrs(seast) = GridElem(irev-1,1,5)%number
1870 44544 : GridElem(i,1,3)%nbrs_face(seast) = 5
1871 44544 : GridElem(i,1,3)%nbrs_wgt(seast) = CornerWgt
1872 44544 : GridElem(i,1,5)%nbrs(seast) = GridElem(irev-1,1,3)%number
1873 44544 : GridElem(i,1,5)%nbrs_face(seast) = 3
1874 44544 : GridElem(i,1,5)%nbrs_wgt(seast) = CornerWgt
1875 : endif
1876 : end do
1877 :
1878 : ! ==================================
1879 : ! south edge of 4 / west edge of 5
1880 : ! ==================================
1881 :
1882 47616 : do i=1,ne
1883 46080 : irev=ne+1-i
1884 46080 : nbrs_used(i,1,4,south) = .true.
1885 46080 : nbrs_used(1,i,5,west) = .true.
1886 :
1887 46080 : GridElem(i,1,4)%nbrs(south) = GridElem(1,i,5)%number
1888 46080 : GridElem(i,1,4)%nbrs_face(south) = 5
1889 46080 : GridElem(i,1,4)%nbrs_wgt(south) = EdgeWgtP
1890 46080 : GridElem(1,i,5)%nbrs(west) = GridElem(i,1,4)%number
1891 46080 : GridElem(1,i,5)%nbrs_face(west) = 4
1892 46080 : GridElem(1,i,5)%nbrs_wgt(west) = EdgeWgtP
1893 : ! Special rules for corner 'edges'
1894 46080 : if( i /= 1) then
1895 44544 : nbrs_used(i,1,4,swest) = .true.
1896 44544 : nbrs_used(1,i,5,swest) = .true.
1897 :
1898 44544 : GridElem(i,1,4)%nbrs(swest) = GridElem(1,i-1,5)%number
1899 44544 : GridElem(i,1,4)%nbrs_face(swest) = 5
1900 44544 : GridElem(i,1,4)%nbrs_wgt(swest) = CornerWgt
1901 44544 : GridElem(1,i,5)%nbrs(swest) = GridElem(i-1,1,4)%number
1902 44544 : GridElem(1,i,5)%nbrs_face(swest) = 4
1903 44544 : GridElem(1,i,5)%nbrs_wgt(swest) = CornerWgt
1904 : endif
1905 47616 : if( i /= ne) then
1906 44544 : nbrs_used(i,1,4,seast) = .true.
1907 44544 : nbrs_used(1,i,5,nwest) = .true.
1908 :
1909 44544 : GridElem(i,1,4)%nbrs(seast) = GridElem(1,i+1,5)%number
1910 44544 : GridElem(i,1,4)%nbrs_face(seast) = 5
1911 44544 : GridElem(i,1,4)%nbrs_wgt(seast) = CornerWgt
1912 44544 : GridElem(1,i,5)%nbrs(nwest) = GridElem(i+1,1,4)%number
1913 44544 : GridElem(1,i,5)%nbrs_face(nwest) = 4
1914 44544 : GridElem(1,i,5)%nbrs_wgt(nwest) = CornerWgt
1915 : endif
1916 : end do
1917 :
1918 : ! ==================================
1919 : ! north edge of 1 / south edge of 6
1920 : ! ==================================
1921 :
1922 47616 : do i=1,ne
1923 46080 : nbrs_used(i,ne,1,north) = .true.
1924 46080 : nbrs_used(i,1,6,south) = .true.
1925 :
1926 :
1927 46080 : GridElem(i,ne,1)%nbrs(north) = GridElem(i,1 ,6)%number
1928 46080 : GridElem(i,ne,1)%nbrs_face(north) = 6
1929 46080 : GridElem(i,ne,1)%nbrs_wgt(north) = EdgeWgtP
1930 46080 : GridElem(i,1 ,6)%nbrs(south) = GridElem(i,ne,1)%number
1931 46080 : GridElem(i,1 ,6)%nbrs_face(south) = 1
1932 46080 : GridElem(i,1 ,6)%nbrs_wgt(south) = EdgeWgtP
1933 : ! Special rules for corner 'edges'
1934 46080 : if( i /= 1) then
1935 44544 : nbrs_used(i,ne,1,nwest) = .true.
1936 44544 : nbrs_used(i,1,6,swest) = .true.
1937 :
1938 44544 : GridElem(i,ne,1)%nbrs(nwest) = GridElem(i-1,1 ,6)%number
1939 44544 : GridElem(i,ne,1)%nbrs_face(nwest) = 6
1940 44544 : GridElem(i,ne,1)%nbrs_wgt(nwest) = CornerWgt
1941 44544 : GridElem(i,1 ,6)%nbrs(swest) = GridElem(i-1,ne,1)%number
1942 44544 : GridElem(i,1 ,6)%nbrs_face(swest) = 1
1943 44544 : GridElem(i,1 ,6)%nbrs_wgt(swest) = CornerWgt
1944 : endif
1945 47616 : if( i /= ne) then
1946 44544 : nbrs_used(i,ne,1,neast) = .true.
1947 44544 : nbrs_used(i,1,6,seast) = .true.
1948 :
1949 :
1950 44544 : GridElem(i,ne,1)%nbrs(neast) = GridElem(i+1,1 ,6)%number
1951 44544 : GridElem(i,ne,1)%nbrs_face(neast) = 6
1952 44544 : GridElem(i,ne,1)%nbrs_wgt(neast) = CornerWgt
1953 44544 : GridElem(i,1 ,6)%nbrs(seast) = GridElem(i+1,ne,1)%number
1954 44544 : GridElem(i,1 ,6)%nbrs_face(seast) = 1
1955 44544 : GridElem(i,1 ,6)%nbrs_wgt(seast) = CornerWgt
1956 : endif
1957 : end do
1958 :
1959 : ! ==================================
1960 : ! north edge of 2 / east edge of 6
1961 : ! ==================================
1962 :
1963 47616 : do i=1,ne
1964 46080 : nbrs_used(i,ne,2,north) = .true.
1965 46080 : nbrs_used(ne,i,6,east ) = .true.
1966 :
1967 46080 : GridElem(i,ne,2)%nbrs(north) = GridElem(ne,i,6)%number
1968 46080 : GridElem(i,ne,2)%nbrs_face(north) = 6
1969 46080 : GridElem(i,ne,2)%nbrs_wgt(north) = EdgeWgtP
1970 46080 : GridElem(ne,i,6)%nbrs(east) = GridElem(i,ne,2)%number
1971 46080 : GridElem(ne,i,6)%nbrs_face(east) = 2
1972 46080 : GridElem(ne,i,6)%nbrs_wgt(east) = EdgeWgtP
1973 : ! Special rules for corner 'edges'
1974 46080 : if( i /= 1) then
1975 44544 : nbrs_used(i,ne,2,nwest) = .true.
1976 44544 : nbrs_used(ne,i,6,seast) = .true.
1977 :
1978 44544 : GridElem(i,ne,2)%nbrs(nwest) = GridElem(ne,i-1,6)%number
1979 44544 : GridElem(i,ne,2)%nbrs_face(nwest) = 6
1980 44544 : GridElem(i,ne,2)%nbrs_wgt(nwest) = CornerWgt
1981 44544 : GridElem(ne,i,6)%nbrs(seast) = GridElem(i-1,ne,2)%number
1982 44544 : GridElem(ne,i,6)%nbrs_face(seast) = 2
1983 44544 : GridElem(ne,i,6)%nbrs_wgt(seast) = CornerWgt
1984 : endif
1985 47616 : if( i /= ne) then
1986 44544 : nbrs_used(i,ne,2,neast) = .true.
1987 44544 : nbrs_used(ne,i,6,neast) = .true.
1988 :
1989 :
1990 44544 : GridElem(i,ne,2)%nbrs(neast) = GridElem(ne,i+1,6)%number
1991 44544 : GridElem(i,ne,2)%nbrs_face(neast) = 6
1992 44544 : GridElem(i,ne,2)%nbrs_wgt(neast) = CornerWgt
1993 44544 : GridElem(ne,i,6)%nbrs(neast) = GridElem(i+1,ne,2)%number
1994 44544 : GridElem(ne,i,6)%nbrs_face(neast) = 2
1995 44544 : GridElem(ne,i,6)%nbrs_wgt(neast) = CornerWgt
1996 : endif
1997 : end do
1998 :
1999 : ! ===================================
2000 : ! north edge of 3 / north edge of 6
2001 : ! ===================================
2002 :
2003 47616 : do i=1,ne
2004 46080 : irev=ne+1-i
2005 46080 : nbrs_used(i,ne,3,north) = .true.
2006 46080 : nbrs_used(i,ne,6,north) = .true.
2007 :
2008 46080 : GridElem(i,ne,3)%nbrs(north) = GridElem(irev,ne,6)%number
2009 46080 : GridElem(i,ne,3)%nbrs_face(north) = 6
2010 46080 : GridElem(i,ne,3)%nbrs_wgt(north) = EdgeWgtP
2011 46080 : GridElem(i,ne,6)%nbrs(north) = GridElem(irev,ne,3)%number
2012 46080 : GridElem(i,ne,6)%nbrs_face(north) = 3
2013 46080 : GridElem(i,ne,6)%nbrs_wgt(north) = EdgeWgtP
2014 : ! Special rules for corner 'edges'
2015 46080 : if( i /= 1) then
2016 44544 : nbrs_used(i,ne,3,nwest) = .true.
2017 44544 : nbrs_used(i,ne,6,nwest) = .true.
2018 :
2019 44544 : GridElem(i,ne,3)%nbrs(nwest) = GridElem(irev+1,ne,6)%number
2020 44544 : GridElem(i,ne,3)%nbrs_face(nwest) = 6
2021 44544 : GridElem(i,ne,3)%nbrs_wgt(nwest) = CornerWgt
2022 44544 : GridElem(i,ne,6)%nbrs(nwest) = GridElem(irev+1,ne,3)%number
2023 44544 : GridElem(i,ne,6)%nbrs_face(nwest) = 3
2024 44544 : GridElem(i,ne,6)%nbrs_wgt(nwest) = CornerWgt
2025 : endif
2026 47616 : if( i /= ne) then
2027 44544 : nbrs_used(i,ne,3,neast) = .true.
2028 44544 : nbrs_used(i,ne,6,neast) = .true.
2029 :
2030 44544 : GridElem(i,ne,3)%nbrs(neast) = GridElem(irev-1,ne,6)%number
2031 44544 : GridElem(i,ne,3)%nbrs_face(neast) = 6
2032 44544 : GridElem(i,ne,3)%nbrs_wgt(neast) = CornerWgt
2033 44544 : GridElem(i,ne,6)%nbrs(neast) = GridElem(irev-1,ne,3)%number
2034 44544 : GridElem(i,ne,6)%nbrs_face(neast) = 3
2035 44544 : GridElem(i,ne,6)%nbrs_wgt(neast) = CornerWgt
2036 : endif
2037 : end do
2038 :
2039 : ! ===================================
2040 : ! north edge of 4 / west edge of 6
2041 : ! ===================================
2042 :
2043 47616 : do i=1,ne
2044 46080 : irev=ne+1-i
2045 46080 : nbrs_used(i,ne,4,north) = .true.
2046 46080 : nbrs_used(1,i,6,west) = .true.
2047 :
2048 46080 : GridElem(i,ne,4)%nbrs(north) = GridElem(1,irev,6)%number
2049 46080 : GridElem(i,ne,4)%nbrs_face(north) = 6
2050 46080 : GridElem(i,ne,4)%nbrs_wgt(north) = EdgeWgtP
2051 46080 : GridElem(1,i,6)%nbrs(west) = GridElem(irev,ne,4)%number
2052 46080 : GridElem(1,i,6)%nbrs_face(west) = 4
2053 46080 : GridElem(1,i,6)%nbrs_wgt(west) = EdgeWgtP
2054 : ! Special rules for corner 'edges'
2055 46080 : if( i /= 1) then
2056 44544 : nbrs_used(i,ne,4,nwest) = .true.
2057 44544 : nbrs_used(1,i,6,swest) = .true.
2058 :
2059 44544 : GridElem(i,ne,4)%nbrs(nwest) = GridElem(1,irev+1,6)%number
2060 44544 : GridElem(i,ne,4)%nbrs_face(nwest) = 6
2061 44544 : GridElem(i,ne,4)%nbrs_wgt(nwest) = CornerWgt
2062 44544 : GridElem(1,i,6)%nbrs(swest) = GridElem(irev+1,ne,4)%number
2063 44544 : GridElem(1,i,6)%nbrs_face(swest) = 4
2064 44544 : GridElem(1,i,6)%nbrs_wgt(swest) = CornerWgt
2065 : endif
2066 47616 : if( i /= ne) then
2067 44544 : nbrs_used(i,ne,4,neast) = .true.
2068 44544 : nbrs_used(1,i,6,nwest) = .true.
2069 :
2070 44544 : GridElem(i,ne,4)%nbrs(neast) = GridElem(1,irev-1,6)%number
2071 44544 : GridElem(i,ne,4)%nbrs_face(neast) = 6
2072 44544 : GridElem(i,ne,4)%nbrs_wgt(neast) = CornerWgt
2073 44544 : GridElem(1,i,6)%nbrs(nwest) = GridElem(irev-1,ne,4)%number
2074 44544 : GridElem(1,i,6)%nbrs_face(nwest) = 4
2075 44544 : GridElem(1,i,6)%nbrs_wgt(nwest) = CornerWgt
2076 : endif
2077 : end do
2078 :
2079 :
2080 : ielem = 1 ! Element counter
2081 10752 : do k=1,6
2082 287232 : do j=1,ne
2083 8580096 : do i=1,ne
2084 8294400 : GridVertex(ielem)%nbrs_ptr(1) = 1
2085 74649600 : do ll=1,8
2086 66355200 : loc = GridVertex(ielem)%nbrs_ptr(ll)
2087 74649600 : if (nbrs_used(i,j,k,ll)) then
2088 66318336 : GridVertex(ielem)%nbrs(loc) = GridElem(i,j,k)%nbrs(ll)
2089 66318336 : GridVertex(ielem)%nbrs_face(loc) = GridElem(i,j,k)%nbrs_face(ll)
2090 66318336 : GridVertex(ielem)%nbrs_wgt(loc) = GridElem(i,j,k)%nbrs_wgt(ll)
2091 66318336 : GridVertex(ielem)%nbrs_wgt_ghost(loc) = GridElem(i,j,k)%nbrs_wgt_ghost(ll)
2092 :
2093 66318336 : GridVertex(ielem)%nbrs_ptr(ll+1) = GridVertex(ielem)%nbrs_ptr(ll)+1
2094 : else
2095 36864 : GridVertex(ielem)%nbrs_ptr(ll+1) = GridVertex(ielem)%nbrs_ptr(ll)
2096 : end if
2097 : end do
2098 8294400 : GridVertex(ielem)%number = GridElem(i,j,k)%number
2099 8294400 : GridVertex(ielem)%processor_number = 0
2100 8294400 : GridVertex(ielem)%SpaceCurve = GridElem(i,j,k)%SpaceCurve
2101 8570880 : ielem=ielem+1
2102 : end do
2103 : end do
2104 : end do
2105 :
2106 1536 : DEALLOCATE(Mesh)
2107 10752 : do k = 1, nfaces
2108 287232 : do j = 1, ne
2109 8580096 : do i = 1, ne
2110 8570880 : call deallocate_gridvertex_nbrs(GridElem(i,j,k))
2111 : end do
2112 : end do
2113 : end do
2114 1536 : DEALLOCATE(GridElem)
2115 1536 : DEALLOCATE(nbrs_used)
2116 :
2117 : ! =======================================
2118 : ! Generate cube graph...
2119 : ! =======================================
2120 :
2121 : ! ============================================
2122 : ! Setup the Grid edges (topology independent)
2123 : ! ============================================
2124 1536 : call initgridedge(GridEdge,GridVertex)
2125 :
2126 : ! ============================================
2127 : ! Setup the Grid edge Indirect addresses
2128 : ! (topology dependent)
2129 : ! ============================================
2130 1536 : nedge = SIZE(GridEdge)
2131 66319872 : do i=1,nedge
2132 66319872 : call CubeSetupEdgeIndex(GridEdge(i))
2133 : enddo
2134 :
2135 1536 : end subroutine CubeTopology
2136 :
2137 : ! ===================================================================
2138 : ! CubeEdgeCount:
2139 : !
2140 : ! Determine the number of Grid Edges
2141 : !
2142 : ! ===================================================================
2143 :
2144 1536 : function CubeEdgeCount() result(nedge)
2145 : use dimensions_mod, only : ne
2146 : implicit none
2147 : integer :: nedge
2148 :
2149 1536 : if (0==ne) call endrun('Error in CubeEdgeCount: ne is zero')
2150 1536 : nedge = nfaces*(ne*ne*nInnerElemEdge - nCornerElemEdge)
2151 :
2152 1536 : end function CubeEdgeCount
2153 :
2154 : ! ===================================================================
2155 : ! CubeElemCount:
2156 : !
2157 : ! Determine the number of Grid Elem
2158 : !
2159 : ! ===================================================================
2160 :
2161 1536 : function CubeElemCount() result(nelem)
2162 :
2163 : use dimensions_mod, only : ne
2164 :
2165 : implicit none
2166 : integer :: nelem
2167 1536 : if (0==ne) call endrun('Error in CubeElemCount: ne is zero')
2168 :
2169 1536 : nelem = nfaces*ne*ne
2170 1536 : end function CubeElemCount
2171 :
2172 66318336 : subroutine CubeSetupEdgeIndex(Edge)
2173 : use gridgraph_mod, only : gridedge_t
2174 : use dimensions_mod, only : np
2175 : use control_mod, only : north, south, east, west, neast, seast, swest, nwest
2176 : type (GridEdge_t),target :: Edge
2177 :
2178 : integer :: np0,sFace,dFace
2179 : logical :: reverse
2180 : integer,allocatable :: forwardV(:), forwardP(:)
2181 : integer,allocatable :: backwardV(:), backwardP(:)
2182 :
2183 66318336 : sFace = Edge%tail_face
2184 66318336 : dFace = Edge%head_face
2185 : ! Do not reverse the indices
2186 66318336 : reverse=.FALSE.
2187 :
2188 : ! Under special conditions use index reversal
2189 : if( (SFace == south .AND. dFace == east) &
2190 : .OR. (sFace == east .AND. dFace == south) &
2191 : .OR. (sFace == north .AND. dFace == west) &
2192 : .OR. (sFace == west .AND. dFace == north) &
2193 : .OR. (sFace == south .AND. dFace == south) &
2194 : .OR. (sFace == north .AND. dFace == north) &
2195 : .OR. (sFace == east .AND. dFace == east ) &
2196 66318336 : .OR. (sFace == west .AND. dFace == west ) ) then
2197 368640 : reverse=.TRUE.
2198 368640 : Edge%reverse=.TRUE.
2199 : endif
2200 :
2201 :
2202 66318336 : end subroutine CubeSetupEdgeIndex
2203 :
2204 : !
2205 : ! HOMME mapping from sphere (or other manifold) to reference element
2206 : ! one should be able to add any mapping here. For each new map,
2207 : ! an associated dmap() routine (which computes the map derivative matrix)
2208 : ! must also be written
2209 : ! Note that for conservation, the parameterization of element edges must be
2210 : ! identical for adjacent elements. (this is violated with HOMME's default
2211 : ! equi-angular cubed-sphere mapping for non-cubed sphere grids, hence the
2212 : ! need for a new map)
2213 : !
2214 345600 : function ref2sphere(a,b, corners3D, ref_map, corners, facenum) result(sphere)
2215 : real(kind=r8) :: a,b
2216 : type (spherical_polar_t) :: sphere
2217 : type (cartesian3d_t) :: corners3D(4)
2218 : integer :: ref_map
2219 : ! only needed for gnominic maps
2220 : type (cartesian2d_t), optional :: corners(4)
2221 : integer, optional :: facenum
2222 :
2223 :
2224 345600 : if (ref_map==0) then
2225 345600 : if (.not. present(corners) ) &
2226 0 : call endrun('ref2sphere(): missing arguments for equiangular map')
2227 345600 : sphere = ref2sphere_equiangular(a,b,corners,facenum)
2228 0 : elseif (ref_map==1) then
2229 0 : call endrun('gnomonic map not yet coded')
2230 0 : elseif (ref_map==2) then
2231 0 : sphere = ref2sphere_elementlocal(a,b,corners3D)
2232 : else
2233 0 : call endrun('ref2sphere(): bad value of ref_map')
2234 : endif
2235 345600 : end function ref2sphere
2236 :
2237 : !
2238 : ! map a point in the referece element to the sphere
2239 : !
2240 345600 : function ref2sphere_equiangular(a,b, corners, face_no) result(sphere)
2241 : implicit none
2242 : real(kind=r8) :: a,b
2243 : integer,intent(in) :: face_no
2244 : type (spherical_polar_t) :: sphere
2245 : type (cartesian2d_t) :: corners(4)
2246 : ! local
2247 : real(kind=r8) :: pi,pj,qi,qj
2248 : type (cartesian2d_t) :: cart
2249 :
2250 : ! map (a,b) to the [-pi/2,pi/2] equi angular cube face: x1,x2
2251 : ! a = gp%points(i)
2252 : ! b = gp%points(j)
2253 345600 : pi = (1-a)/2
2254 345600 : pj = (1-b)/2
2255 345600 : qi = (1+a)/2
2256 345600 : qj = (1+b)/2
2257 : cart%x = pi*pj*corners(1)%x &
2258 : + qi*pj*corners(2)%x &
2259 : + qi*qj*corners(3)%x &
2260 345600 : + pi*qj*corners(4)%x
2261 : cart%y = pi*pj*corners(1)%y &
2262 : + qi*pj*corners(2)%y &
2263 : + qi*qj*corners(3)%y &
2264 345600 : + pi*qj*corners(4)%y
2265 : ! map from [pi/2,pi/2] equ angular cube face to sphere:
2266 345600 : sphere=projectpoint(cart,face_no)
2267 :
2268 345600 : end function ref2sphere_equiangular
2269 :
2270 : !-----------------------------------------------------------------------------------------
2271 : ! ELEMENT LOCAL MAP (DOES NOT USE CUBE FACES)
2272 : ! unlike gnomonic equiangular map, this map will map all straight lines to
2273 : ! great circle arcs
2274 : !
2275 : ! map a point in the referece element to the quad on the sphere by a
2276 : ! general map, without using faces the map works this way: first, fix
2277 : ! a coordinate (say, X). Map 4 corners of the ref element (corners are
2278 : ! (-1,-1),(-1,1),(1,1), and (1,-1)) into 4 X-components of the quad in
2279 : ! physical space via a bilinear map. Do so for Y and Z components as
2280 : ! well. It produces a map: Ref element (\xi, \eta) ---> A quad in XYZ
2281 : ! (ess, a piece of a twisted plane) with vertices of our target quad. though
2282 : ! the quad lies in a plane and not on the sphere manifold, its
2283 : ! vertices belong to the sphere (by initial conditions). The last step
2284 : ! is to utilize a map (X,Y,X) --> (X,Y,Z)/SQRT(X**2+Y**2+Z**2) to
2285 : ! project the quad to the unit sphere.
2286 : ! -----------------------------------------------------------------------------------------
2287 0 : function ref2sphere_elementlocal(a,b, corners3D) result(sphere)
2288 : use element_mod, only : element_t
2289 : implicit none
2290 : real(kind=r8) :: a,b
2291 : type (cartesian3d_t) :: corners3D(4)
2292 : type (spherical_polar_t) :: sphere
2293 : real(kind=r8) :: q(4) ! local
2294 :
2295 0 : q(1)=(1-a)*(1-b); q(2)=(1+a)*(1-b); q(3)=(1+a)*(1+b); q(4)=(1-a)*(1+b);
2296 0 : q=q/4.0_r8;
2297 0 : sphere = ref2sphere_elementlocal_q(q,corners3D)
2298 0 : end function ref2sphere_elementlocal
2299 :
2300 0 : function ref2sphere_elementlocal_q(q, corners) result(sphere)
2301 : implicit none
2302 : real(kind=r8) :: q(4)
2303 : type (spherical_polar_t) :: sphere
2304 : type (cartesian3d_t) :: corners(4)
2305 : ! local
2306 : type (cartesian3d_t) :: cart
2307 : real(kind=r8) :: c(3,4), xx(3), r
2308 : integer :: i
2309 :
2310 : !3D corners fo the quad
2311 0 : c(1,1)=corners(1)%x; c(2,1)=corners(1)%y; c(3,1)=corners(1)%z;
2312 0 : c(1,2)=corners(2)%x; c(2,2)=corners(2)%y; c(3,2)=corners(2)%z;
2313 0 : c(1,3)=corners(3)%x; c(2,3)=corners(3)%y; c(3,3)=corners(3)%z;
2314 0 : c(1,4)=corners(4)%x; c(2,4)=corners(4)%y; c(3,4)=corners(4)%z;
2315 :
2316 : !physical point on a plane (sliced), not yet on the sphere
2317 0 : do i=1,3
2318 0 : xx(i)=sum(c(i,:)*q(:))
2319 : end do
2320 :
2321 : !distance from the plane point to the origin
2322 0 : r = sqrt(xx(1)**2+xx(2)**2+xx(3)**2)
2323 :
2324 : !projecting the plane point to the sphere
2325 0 : cart%x=xx(1)/r; cart%y=xx(2)/r; cart%z=xx(3)/r;
2326 :
2327 : !XYZ coords of the point to lon/lat
2328 0 : sphere=change_coordinates(cart)
2329 :
2330 0 : end function ref2sphere_elementlocal_q
2331 :
2332 0 : end module cube_mod
|
