Line data Source code
1 : module element_mod
2 :
3 : use shr_kind_mod, only: r8=>shr_kind_r8, i8=>shr_kind_i8
4 : use coordinate_systems_mod, only: spherical_polar_t, cartesian2D_t, cartesian3D_t, distance
5 : use dimensions_mod, only: np, nc, npsq, nlev, nlevp, qsize_d, max_neigh_edges,ntrac_d
6 : use edgetype_mod, only: edgedescriptor_t
7 : use gridgraph_mod, only: gridvertex_t
8 :
9 : implicit none
10 : private
11 : integer, public, parameter :: timelevels = 3
12 :
13 :
14 : ! =========== PRIMITIVE-EQUATION DATA-STRUCTURES =====================
15 :
16 : type, public :: elem_state_t
17 :
18 : ! prognostic variables for preqx solver
19 :
20 : ! prognostics must match those in prim_restart_mod.F90
21 : ! vertically-lagrangian code advects dp3d instead of ps
22 : ! tracers Q, Qdp always use 2 level time scheme
23 :
24 : real (kind=r8) :: v (np,np,2,nlev,timelevels) ! velocity
25 : real (kind=r8) :: T (np,np,nlev,timelevels) ! temperature
26 : real (kind=r8) :: dp3d (np,np,nlev,timelevels) ! dry delta p on levels
27 : real (kind=r8) :: psdry (np,np) ! dry surface pressure
28 : real (kind=r8) :: phis (np,np) ! surface geopotential (prescribed)
29 : real (kind=r8), allocatable :: Qdp(:,:,:,:,:) ! Tracer mass
30 : end type elem_state_t
31 :
32 : !___________________________________________________________________
33 : type, public :: derived_state_t
34 : !
35 : ! storage for subcycling tracers/dynamics
36 : !
37 : real (kind=r8) :: vn0 (np,np,2,nlev) ! velocity for SE tracer advection
38 : real (kind=r8) :: dpdiss_biharmonic(np,np,nlev) ! mean dp dissipation tendency, if nu_p>0
39 : real (kind=r8) :: dpdiss_ave(np,np,nlev) ! mean dp used to compute psdiss_tens
40 :
41 : ! diagnostics for explicit timestep
42 : real (kind=r8) :: phi(np,np,nlev) ! geopotential
43 : real (kind=r8) :: omega(np,np,nlev) ! vertical velocity
44 :
45 : ! tracer advection fields used for consistency and limiters
46 : real (kind=r8) :: dp(np,np,nlev) ! for dp_tracers at physics timestep
47 : real (kind=r8), allocatable :: divdp(:,:,:) ! divergence of dp
48 : real (kind=r8), allocatable :: divdp_proj(:,:,:) ! DSSed divdp
49 : real (kind=r8) :: mass(MAX(qsize_d,ntrac_d)+9) ! total tracer mass for diagnostics
50 :
51 : ! forcing terms for CAM
52 : real (kind=r8), allocatable :: FQ(:,:,:,:) ! tracer forcing
53 : real (kind=r8) :: FM(np,np,2,nlev) ! momentum forcing
54 : real (kind=r8), allocatable :: FDP(:,:,:) ! save full updated dp right after physics
55 : real (kind=r8) :: FT(np,np,nlev) ! temperature forcing
56 : real (kind=r8) :: etadot_prescribed(np,np,nlevp) ! prescribed vertical tendency
57 : real (kind=r8) :: u_met(np,np,nlev) ! zonal component of prescribed meteorology winds
58 : real (kind=r8) :: dudt_met(np,np,nlev) ! rate of change of zonal component of prescribed meteorology winds
59 : real (kind=r8) :: v_met(np,np,nlev) ! meridional component of prescribed meteorology winds
60 : real (kind=r8) :: dvdt_met(np,np,nlev) ! rate of change of meridional component of prescribed meteorology winds
61 : real (kind=r8) :: T_met(np,np,nlev) ! prescribed meteorology temperature
62 : real (kind=r8) :: dTdt_met(np,np,nlev) ! rate of change of prescribed meteorology temperature
63 : real (kind=r8) :: ps_met(np,np) ! surface pressure of prescribed meteorology
64 : real (kind=r8) :: dpsdt_met(np,np) ! rate of change of surface pressure of prescribed meteorology
65 : real (kind=r8) :: nudge_factor(np,np,nlev) ! nudging factor (prescribed)
66 : real (kind=r8) :: Utnd(npsq,nlev) ! accumulated U tendency due to nudging towards prescribed met
67 : real (kind=r8) :: Vtnd(npsq,nlev) ! accumulated V tendency due to nudging towards prescribed met
68 : real (kind=r8) :: Ttnd(npsq,nlev) ! accumulated T tendency due to nudging towards prescribed met
69 :
70 : real (kind=r8) :: pecnd(np,np,nlev) ! pressure perturbation from condensate
71 :
72 : ! reference profiles
73 : real (kind=r8) :: T_ref(np,np,nlev) ! reference temperature
74 : real (kind=r8) :: dp_ref(np,np,nlev) ! reference pressure level thickness
75 : end type derived_state_t
76 :
77 : !___________________________________________________________________
78 : type, public :: elem_accum_t
79 :
80 :
81 : ! the "4" timelevels represents data computed at:
82 : ! 1 t-.5
83 : ! 2 t+.5 after dynamics
84 : ! 3 t+.5 after forcing
85 : ! 4 t+.5 after Robert
86 : ! after calling TimeLevelUpdate, all times above decrease by 1.0
87 :
88 :
89 : end type elem_accum_t
90 :
91 :
92 : ! ============= DATA-STRUCTURES COMMON TO ALL SOLVERS ================
93 :
94 : type, public :: index_t
95 : integer :: ia(npsq),ja(npsq)
96 : integer :: is,ie
97 : integer :: NumUniquePts
98 : integer :: UniquePtOffset
99 : end type index_t
100 :
101 : !___________________________________________________________________
102 : type, public :: element_t
103 : integer :: LocalId
104 : integer :: GlobalId
105 :
106 : ! Coordinate values of element points
107 : type (spherical_polar_t) :: spherep(np,np) ! Spherical coords of GLL points
108 :
109 : ! Equ-angular gnomonic projection coordinates
110 : type (cartesian2D_t) :: cartp(np,np) ! gnomonic coords of GLL points
111 : type (cartesian2D_t) :: corners(4) ! gnomonic coords of element corners
112 : real (kind=r8) :: u2qmap(4,2) ! bilinear map from ref element to quad in cubedsphere coordinates
113 : ! SHOULD BE REMOVED
114 : ! 3D cartesian coordinates
115 : type (cartesian3D_t) :: corners3D(4)
116 :
117 : ! Element diagnostics
118 : real (kind=r8) :: area ! Area of element
119 : real (kind=r8) :: normDinv ! some type of norm of Dinv used for CFL
120 : real (kind=r8) :: dx_short ! short length scale in km
121 : real (kind=r8) :: dx_long ! long length scale in km
122 :
123 : real (kind=r8) :: variable_hyperviscosity(np,np) ! hyperviscosity based on above
124 : real (kind=r8) :: hv_courant ! hyperviscosity courant number
125 : real (kind=r8) :: tensorVisc(np,np,2,2) !og, matrix V for tensor viscosity
126 :
127 : ! Edge connectivity information
128 : ! integer :: node_numbers(4)
129 : ! integer :: node_multiplicity(4) ! number of elements sharing corner node
130 :
131 : type (GridVertex_t) :: vertex ! element grid vertex information
132 : type (EdgeDescriptor_t) :: desc
133 :
134 : type (elem_state_t) :: state
135 :
136 : type (derived_state_t) :: derived
137 : ! Metric terms
138 : real (kind=r8) :: met(np,np,2,2) ! metric tensor on velocity and pressure grid
139 : real (kind=r8) :: metinv(np,np,2,2) ! metric tensor on velocity and pressure grid
140 : real (kind=r8) :: metdet(np,np) ! g = SQRT(det(g_ij)) on velocity and pressure grid
141 : real (kind=r8) :: rmetdet(np,np) ! 1/metdet on velocity pressure grid
142 : real (kind=r8) :: D(np,np,2,2) ! Map covariant field on cube to vector field on the sphere
143 : real (kind=r8) :: Dinv(np,np,2,2) ! Map vector field on the sphere to covariant v on cube
144 :
145 :
146 : ! Mass flux across the sides of each sub-element.
147 : ! The storage is redundent since the mass across shared sides
148 : ! must be equal in magnitude and opposite in sign.
149 : ! The layout is like:
150 : ! --------------------------------------------------------------
151 : ! ^| (1,4,3) | | | (4,4,3) |
152 : ! || | | | |
153 : ! ||(1,4,4) | | |(4,4,4) |
154 : ! || (1,4,2)| | | (4,4,2)|
155 : ! || | | | |
156 : ! || (1,4,1) | | | (4,4,1) |
157 : ! |---------------------------------------------------------------
158 : ! S| | | | |
159 : ! e| | | | |
160 : ! c| | | | |
161 : ! o| | | | |
162 : ! n| | | | |
163 : ! d| | | | |
164 : ! ---------------------------------------------------------------
165 : ! C| | | | |
166 : ! o| | | | |
167 : ! o| | | | |
168 : ! r| | | | |
169 : ! d| | | | |
170 : ! i| | | | |
171 : ! n---------------------------------------------------------------
172 : ! a| (1,1,3) | | | (4,1,3) |
173 : ! t| | | |(4,1,4) |
174 : ! e|(1,1,4) | | | |
175 : ! | (1,1,2)| | | (4,1,2)|
176 : ! | | | | |
177 : ! | (1,1,1) | | | (4,1,1) |
178 : ! ---------------------------------------------------------------
179 : ! First Coordinate ------->
180 : real (kind=r8) :: sub_elem_mass_flux(nc,nc,4,nlev)
181 :
182 : ! Convert vector fields from spherical to rectangular components
183 : ! The transpose of this operation is its pseudoinverse.
184 : real (kind=r8) :: vec_sphere2cart(np,np,3,2)
185 :
186 : ! Mass matrix terms for an element on a cube face
187 : real (kind=r8) :: mp(np,np) ! mass matrix on v and p grid
188 : real (kind=r8) :: rmp(np,np) ! inverse mass matrix on v and p grid
189 :
190 : ! Mass matrix terms for an element on the sphere
191 : ! This mass matrix is used when solving the equations in weak form
192 : ! with the natural (surface area of the sphere) inner product
193 : real (kind=r8) :: spheremp(np,np) ! mass matrix on v and p grid
194 : real (kind=r8) :: rspheremp(np,np) ! inverse mass matrix on v and p grid
195 :
196 : integer(i8) :: gdofP(np,np) ! global degree of freedom (P-grid)
197 :
198 : real (kind=r8) :: fcor(np,np) ! Coreolis term
199 :
200 : type (index_t) :: idxP
201 : type (index_t),pointer :: idxV
202 : integer :: FaceNum
203 :
204 : ! force element_t to be a multiple of 8 bytes.
205 : ! on BGP, code will crash (signal 7, or signal 15) if 8 byte alignment is off
206 : ! check core file for:
207 : ! core.63:Generated by interrupt..(Alignment Exception DEAR=0xa1ef671c ESR=0x01800000 CCR0=0x4800a002)
208 : integer :: dummy
209 : end type element_t
210 :
211 : !___________________________________________________________________
212 : public :: element_coordinates
213 : public :: element_var_coordinates
214 : public :: element_var_coordinates3D
215 : public :: GetColumnIdP,GetColumnIdV
216 : public :: allocate_element_desc
217 : public :: PrintElem
218 :
219 : contains
220 :
221 0 : subroutine PrintElem(arr)
222 :
223 : real(kind=r8) :: arr(:,:)
224 : integer :: i,j
225 :
226 0 : do j=np,1,-1
227 0 : write(6,*) (arr(i,j), i=1,np)
228 : enddo
229 :
230 0 : end subroutine PrintElem
231 : ! ===================== ELEMENT_MOD METHODS ==========================
232 :
233 0 : function GetColumnIdP(elem,i,j) result(col_id)
234 :
235 : ! Get unique identifier for a Physics column on the P-grid
236 :
237 : type(element_t), intent(in) :: elem
238 : integer, intent(in) :: i,j
239 : integer :: col_id
240 0 : col_id = elem%gdofP(i,j)
241 0 : end function GetColumnIdP
242 :
243 : !___________________________________________________________________
244 0 : function GetColumnIdV(elem,i,j) result(col_id)
245 :
246 : ! Get unique identifier for a Physics column on the V-grid
247 :
248 : type(element_t), intent(in) :: elem
249 : integer, intent(in) :: i,j
250 : integer :: col_id
251 0 : col_id = elem%gdofP(i,j)
252 0 : end function GetColumnIdV
253 :
254 : !___________________________________________________________________
255 0 : function element_coordinates(start,end,points) result(cart)
256 :
257 : ! Initialize 2D rectilinear element colocation points
258 :
259 : type (cartesian2D_t), intent(in) :: start
260 : type (cartesian2D_t), intent(in) :: end
261 : real(r8), intent(in) :: points(:)
262 : type (cartesian2D_t) :: cart(SIZE(points),SIZE(points))
263 :
264 : type (cartesian2D_t) :: length, centroid
265 : real(r8) :: y
266 : integer :: i,j
267 :
268 0 : length%x = 0.50D0*(end%x-start%x)
269 0 : length%y = 0.50D0*(end%y-start%y)
270 0 : centroid%x = 0.50D0*(end%x+start%x)
271 0 : centroid%y = 0.50D0*(end%y+start%y)
272 0 : do j=1,SIZE(points)
273 0 : y = centroid%y + length%y*points(j)
274 0 : do i=1,SIZE(points)
275 0 : cart(i,j)%x = centroid%x + length%x*points(i)
276 0 : cart(i,j)%y = y
277 : end do
278 : end do
279 0 : end function element_coordinates
280 :
281 : !___________________________________________________________________
282 21600 : function element_var_coordinates(c,points) result(cart)
283 :
284 : type (cartesian2D_t), intent(in) :: c(4)
285 : real(r8), intent(in) :: points(:)
286 : type (cartesian2D_t) :: cart(SIZE(points),SIZE(points))
287 :
288 43200 : real(r8) :: p(size(points))
289 43200 : real(r8) :: q(size(points))
290 : integer :: i,j
291 :
292 129600 : p(:) = (1.0D0-points(:))/2.0D0
293 108000 : q(:) = (1.0D0+points(:))/2.0D0
294 :
295 108000 : do j=1,SIZE(points)
296 453600 : do i=1,SIZE(points)
297 691200 : cart(i,j)%x = p(i)*p(j)*c(1)%x &
298 : + q(i)*p(j)*c(2)%x &
299 : + q(i)*q(j)*c(3)%x &
300 691200 : + p(i)*q(j)*c(4)%x
301 : cart(i,j)%y = p(i)*p(j)*c(1)%y &
302 : + q(i)*p(j)*c(2)%y &
303 : + q(i)*q(j)*c(3)%y &
304 432000 : + p(i)*q(j)*c(4)%y
305 : end do
306 : end do
307 21600 : end function element_var_coordinates
308 :
309 : !___________________________________________________________________
310 0 : function element_var_coordinates3d(c,points) result(cart)
311 :
312 : type(cartesian3D_t), intent(in) :: c(4)
313 : real(r8), intent(in) :: points(:)
314 :
315 : type(cartesian3D_t) :: cart(SIZE(points),SIZE(points))
316 :
317 0 : real(r8) :: p(size(points))
318 0 : real(r8) :: q(size(points)), r
319 : integer :: i,j
320 :
321 0 : p(:) = (1.0D0-points(:))/2.0D0
322 0 : q(:) = (1.0D0+points(:))/2.0D0
323 :
324 0 : do j=1,SIZE(points)
325 0 : do i=1,SIZE(points)
326 0 : cart(i,j)%x = p(i)*p(j)*c(1)%x &
327 : + q(i)*p(j)*c(2)%x &
328 : + q(i)*q(j)*c(3)%x &
329 0 : + p(i)*q(j)*c(4)%x
330 : cart(i,j)%y = p(i)*p(j)*c(1)%y &
331 : + q(i)*p(j)*c(2)%y &
332 : + q(i)*q(j)*c(3)%y &
333 0 : + p(i)*q(j)*c(4)%y
334 : cart(i,j)%z = p(i)*p(j)*c(1)%z &
335 : + q(i)*p(j)*c(2)%z &
336 : + q(i)*q(j)*c(3)%z &
337 0 : + p(i)*q(j)*c(4)%z
338 :
339 : ! project back to sphere:
340 0 : r = distance(cart(i,j))
341 0 : cart(i,j)%x = cart(i,j)%x/r
342 0 : cart(i,j)%y = cart(i,j)%y/r
343 0 : cart(i,j)%z = cart(i,j)%z/r
344 : end do
345 : end do
346 0 : end function element_var_coordinates3d
347 :
348 : !___________________________________________________________________
349 1536 : subroutine allocate_element_desc(elem)
350 :
351 : type (element_t), intent(inout) :: elem(:)
352 : integer :: num, j,i
353 :
354 1536 : num = SIZE(elem)
355 :
356 12336 : do j=1,num
357 32400 : allocate(elem(j)%desc%putmapP(max_neigh_edges))
358 21600 : allocate(elem(j)%desc%getmapP(max_neigh_edges))
359 21600 : allocate(elem(j)%desc%putmapP_ghost(max_neigh_edges))
360 21600 : allocate(elem(j)%desc%getmapP_ghost(max_neigh_edges))
361 21600 : allocate(elem(j)%desc%putmapS(max_neigh_edges))
362 21600 : allocate(elem(j)%desc%getmapS(max_neigh_edges))
363 21600 : allocate(elem(j)%desc%reverse(max_neigh_edges))
364 21600 : allocate(elem(j)%desc%globalID(max_neigh_edges))
365 21600 : allocate(elem(j)%desc%loc2buf(max_neigh_edges))
366 98736 : do i=1,max_neigh_edges
367 86400 : elem(j)%desc%loc2buf(i)=i
368 97200 : elem(j)%desc%globalID(i)=-1
369 : enddo
370 :
371 : end do
372 1536 : end subroutine allocate_element_desc
373 :
374 :
375 0 : end module element_mod
|
