Line data Source code
1 : module prim_init
2 :
3 : use shr_kind_mod, only: r8=>shr_kind_r8
4 : use dimensions_mod, only: nc, use_cslam
5 : use reduction_mod, only: reductionbuffer_ordered_1d_t
6 : use quadrature_mod, only: quadrature_t, gausslobatto
7 :
8 : implicit none
9 : private
10 : save
11 :
12 : public :: prim_init1
13 :
14 : real(r8), public :: fvm_corners(nc+1) ! fvm cell corners on reference element
15 : real(r8), public :: fvm_points(nc) ! fvm cell centers on reference element
16 :
17 : type (quadrature_t), public :: gp ! element GLL points
18 : type (ReductionBuffer_ordered_1d_t) :: red ! reduction buffer (shared)
19 :
20 : contains
21 1536 : subroutine prim_init1(elem, fvm, par, Tl)
22 : use cam_logfile, only: iulog
23 : use shr_sys_mod, only: shr_sys_flush
24 : use thread_mod, only: max_num_threads
25 : use dimensions_mod, only: np, nlev, nelem, nelemd, nelemdmax, qsize_d
26 : use dimensions_mod, only: GlobalUniqueCols, fv_nphys,irecons_tracer
27 : use control_mod, only: topology, partmethod
28 : use element_mod, only: element_t, allocate_element_desc
29 : use fvm_mod, only: fvm_init1
30 : use mesh_mod, only: MeshUseMeshFile
31 : use se_dyn_time_mod, only: timelevel_init, timelevel_t
32 : use mass_matrix_mod, only: mass_matrix
33 : use derivative_mod, only: allocate_subcell_integration_matrix_cslam
34 : use derivative_mod, only: allocate_subcell_integration_matrix_physgrid
35 : use cube_mod, only: cubeedgecount , cubeelemcount, cubetopology
36 : use cube_mod, only: cube_init_atomic, rotation_init_atomic, set_corner_coordinates
37 : use cube_mod, only: assign_node_numbers_to_elem
38 : use mesh_mod, only: MeshSetCoordinates, MeshUseMeshFile, MeshCubeTopology
39 : use mesh_mod, only: MeshCubeElemCount, MeshCubeEdgeCount
40 : use metagraph_mod, only: metavertex_t, localelemcount, initmetagraph, printmetavertex
41 : use gridgraph_mod, only: gridvertex_t, gridedge_t
42 : use gridgraph_mod, only: allocate_gridvertex_nbrs, deallocate_gridvertex_nbrs
43 : use schedtype_mod, only: schedule
44 : use schedule_mod, only: genEdgeSched
45 : use prim_advection_mod, only: prim_advec_init1
46 : use cam_abortutils, only: endrun
47 : use spmd_utils, only: mpi_integer, mpi_max
48 : use parallel_mod, only: parallel_t, syncmp, global_shared_buf, nrepro_vars
49 : use spacecurve_mod, only: genspacepart
50 : use dof_mod, only: global_dof, CreateUniqueIndex, SetElemOffset
51 : use params_mod, only: SFCURVE
52 : use physconst, only: pi
53 : use reduction_mod, only: red_min, red_max, red_max_int, red_flops
54 : use reduction_mod, only: red_sum, red_sum_int, initreductionbuffer
55 : use infnan, only: nan, assignment(=)
56 : use shr_reprosum_mod, only: repro_sum => shr_reprosum_calc
57 : use fvm_analytic_mod, only: compute_basic_coordinate_vars
58 : use fvm_control_volume_mod, only: fvm_struct, allocate_physgrid_vars
59 : use air_composition, only: thermodynamic_active_species_num
60 :
61 : type(element_t), pointer :: elem(:)
62 : type(fvm_struct), pointer :: fvm(:)
63 : type(parallel_t), intent(inout) :: par
64 : type(timelevel_t), intent(out) :: Tl
65 :
66 : ! Local Variables
67 1536 : type (GridVertex_t), target,allocatable :: GridVertex(:)
68 1536 : type (GridEdge_t), target,allocatable :: Gridedge(:)
69 1536 : type (MetaVertex_t), target,allocatable :: MetaVertex(:)
70 :
71 : integer :: ie
72 : integer :: nets, nete
73 : integer :: nelem_edge
74 : integer :: ierr=0, j
75 : logical, parameter :: Debug = .FALSE.
76 :
77 1536 : real(r8), allocatable :: aratio(:,:)
78 : real(r8) :: area(1), xtmp
79 : character(len=80) :: rot_type ! cube edge rotation type
80 :
81 : integer :: i
82 :
83 : character(len=128) :: errmsg
84 : character(len=*), parameter :: subname = 'PRIM_INIT1: '
85 :
86 : ! ====================================
87 : ! Set cube edge rotation type for model
88 : ! unnecessary complication here: all should
89 : ! be on the same footing. RDL
90 : ! =====================================
91 1536 : rot_type = "contravariant"
92 :
93 : ! ===============================================================
94 : ! Allocate and initialize the graph (array of GridVertex_t types)
95 : ! ===============================================================
96 :
97 1536 : if (topology=="cube") then
98 :
99 1536 : if (par%masterproc) then
100 2 : write(iulog,*) subname, "creating cube topology..."
101 2 : call shr_sys_flush(iulog)
102 : end if
103 :
104 1536 : if (MeshUseMeshFile) then
105 0 : nelem = MeshCubeElemCount()
106 0 : nelem_edge = MeshCubeEdgeCount()
107 : else
108 1536 : nelem = CubeElemCount()
109 1536 : nelem_edge = CubeEdgeCount()
110 : end if
111 :
112 8299008 : allocate(GridVertex(nelem))
113 66322944 : allocate(GridEdge(nelem_edge))
114 :
115 8295936 : do j = 1, nelem
116 8295936 : call allocate_gridvertex_nbrs(GridVertex(j))
117 : end do
118 :
119 1536 : if (MeshUseMeshFile) then
120 0 : if (par%masterproc) then
121 0 : write(iulog,*) subname, "Set up grid vertex from mesh..."
122 : end if
123 0 : call MeshCubeTopology(GridEdge, GridVertex)
124 : else
125 1536 : call CubeTopology(GridEdge,GridVertex)
126 : end if
127 :
128 1536 : if (par%masterproc) then
129 2 : write(iulog,*)"...done."
130 : end if
131 : end if
132 1536 : if(par%masterproc) then
133 2 : write(iulog,*) subname, "total number of elements nelem = ",nelem
134 : end if
135 :
136 1536 : if(partmethod == SFCURVE) then
137 1536 : if(par%masterproc) then
138 2 : write(iulog,*) subname, "partitioning graph using SF Curve..."
139 : end if
140 1536 : call genspacepart(GridVertex)
141 : else
142 0 : write(errmsg, *) 'Unsupported partition method, ',partmethod
143 0 : call endrun(subname//trim(errmsg))
144 : end if
145 :
146 : ! ===========================================================
147 : ! given partition, count number of local element descriptors
148 : ! ===========================================================
149 1536 : allocate(MetaVertex(1))
150 1536 : allocate(Schedule(1))
151 :
152 1536 : nelem_edge = SIZE(GridEdge)
153 :
154 : ! ====================================================
155 : ! Generate the communication graph
156 : ! ====================================================
157 1536 : call initMetaGraph(par%rank+1,MetaVertex(1),GridVertex,GridEdge)
158 :
159 1536 : nelemd = LocalElemCount(MetaVertex(1))
160 : if (par%masterproc .and. Debug) then
161 : call PrintMetaVertex(MetaVertex(1))
162 : endif
163 :
164 1536 : if(nelemd <= 0) then
165 0 : call endrun(subname//'Not yet ready to handle nelemd = 0 yet' )
166 : end if
167 1536 : call mpi_allreduce(nelemd, nelemdmax, 1, MPI_INTEGER, MPI_MAX, par%comm, ierr)
168 :
169 : !Allocate elements:
170 1536 : if (nelemd > 0) then
171 15408 : allocate(elem(nelemd))
172 1536 : call allocate_element_desc(elem)
173 : !Allocate Qdp and derived FQ arrays:
174 1536 : if(fv_nphys > 0) then !SE-CSLAM
175 12336 : do ie=1,nelemd
176 32400 : allocate(elem(ie)%state%Qdp(np,np,nlev,thermodynamic_active_species_num,1), stat=ierr)
177 10800 : if( ierr /= 0 ) then
178 0 : call endrun('prim_init1: failed to allocate Qdp array')
179 : end if
180 32400 : allocate(elem(ie)%derived%FQ(np,np,nlev,thermodynamic_active_species_num), stat=ierr)
181 12336 : if( ierr /= 0 ) then
182 0 : call endrun('prim_init1: failed to allocate fq array')
183 : end if
184 : end do
185 : else !Regular SE
186 0 : do ie=1,nelemd
187 0 : allocate(elem(ie)%state%Qdp(np,np,nlev,qsize_d,2), stat=ierr)
188 0 : if( ierr /= 0 ) then
189 0 : call endrun('prim_init1: failed to allocate Qdp array')
190 : end if
191 0 : allocate(elem(ie)%derived%FQ(np,np,nlev,qsize_d), stat=ierr)
192 0 : if( ierr /= 0 ) then
193 0 : call endrun('prim_init1: failed to allocate fq array')
194 : end if
195 : end do
196 : end if
197 : !Allocate remaining derived quantity arrays:
198 12336 : do ie=1,nelemd
199 10800 : allocate(elem(ie)%derived%FDP(np,np,nlev), stat=ierr)
200 10800 : if( ierr /= 0 ) then
201 0 : call endrun('prim_init1: failed to allocate fdp array')
202 : end if
203 10800 : allocate(elem(ie)%derived%divdp(np,np,nlev), stat=ierr)
204 10800 : if( ierr /= 0 ) then
205 0 : call endrun('prim_init1: failed to allocate divdp array')
206 : end if
207 10800 : allocate(elem(ie)%derived%divdp_proj(np,np,nlev), stat=ierr)
208 12336 : if( ierr /= 0 ) then
209 0 : call endrun('prim_init1: failed to allocate divdp_proj array')
210 : end if
211 : end do
212 : end if
213 :
214 1536 : if (fv_nphys > 0) then
215 15408 : allocate(fvm(nelemd))
216 1536 : call allocate_physgrid_vars(fvm,par)
217 : else
218 : ! Even if fvm not needed, still desirable to allocate it as empty
219 : ! so it can be passed as a (size zero) array rather than pointer.
220 0 : allocate(fvm(0))
221 : end if
222 :
223 : ! ====================================================
224 : ! Generate the communication schedule
225 : ! ====================================================
226 :
227 1536 : call genEdgeSched(par, elem, par%rank+1, Schedule(1), MetaVertex(1))
228 :
229 4608 : allocate(global_shared_buf(nelemd, nrepro_vars))
230 12621264 : global_shared_buf = 0.0_r8
231 :
232 1536 : call syncmp(par)
233 :
234 : ! =================================================================
235 : ! Set number of domains (for 'decompose') equal to number of threads
236 : ! for OpenMP across elements, equal to 1 for OpenMP within element
237 : ! =================================================================
238 :
239 : ! =================================================================
240 : ! Initialize shared boundary_exchange and reduction buffers
241 : ! =================================================================
242 1536 : if(par%masterproc) then
243 2 : write(iulog,*) subname, 'init shared boundary_exchange buffers'
244 2 : call shr_sys_flush(iulog)
245 : end if
246 1536 : call InitReductionBuffer(red,3*nlev,max_num_threads)
247 1536 : call InitReductionBuffer(red_sum,5)
248 1536 : call InitReductionBuffer(red_sum_int,1)
249 1536 : call InitReductionBuffer(red_max,1)
250 1536 : call InitReductionBuffer(red_max_int,1)
251 1536 : call InitReductionBuffer(red_min,1)
252 1536 : call initReductionBuffer(red_flops,1)
253 :
254 1536 : gp = gausslobatto(np) ! GLL points
255 :
256 : ! fvm nodes are equally spaced in alpha/beta
257 : ! HOMME with equ-angular gnomonic projection maps alpha/beta space
258 : ! to the reference element via simple scale + translation
259 : ! thus, fvm nodes in reference element [-1,1] are a tensor product of
260 : ! array 'fvm_corners(:)' computed below:
261 1536 : xtmp = nc
262 7680 : do i = 1, nc+1
263 7680 : fvm_corners(i)= 2*(i-1)/xtmp - 1 ! [-1,1] including end points
264 : end do
265 6144 : do i = 1, nc
266 6144 : fvm_points(i)= ( fvm_corners(i)+fvm_corners(i+1) ) /2
267 : end do
268 :
269 1536 : if (topology == "cube") then
270 1536 : if(par%masterproc) then
271 2 : write(iulog,*) subname, "initializing cube elements..."
272 2 : call shr_sys_flush(iulog)
273 : end if
274 1536 : if (MeshUseMeshFile) then
275 0 : call MeshSetCoordinates(elem)
276 : else
277 12336 : do ie = 1, nelemd
278 12336 : call set_corner_coordinates(elem(ie))
279 : end do
280 1536 : call assign_node_numbers_to_elem(elem, GridVertex)
281 : end if
282 12336 : do ie = 1, nelemd
283 12336 : call cube_init_atomic(elem(ie),gp%points)
284 : end do
285 : end if
286 :
287 : ! =================================================================
288 : ! Initialize mass_matrix
289 : ! =================================================================
290 1536 : if(par%masterproc) then
291 2 : write(iulog,*) subname, 'running mass_matrix'
292 2 : call shr_sys_flush(iulog)
293 : end if
294 1536 : call mass_matrix(par, elem)
295 4608 : allocate(aratio(nelemd,1))
296 :
297 1536 : if (topology == "cube") then
298 1536 : area = 0
299 12336 : do ie = 1, nelemd
300 228336 : aratio(ie,1) = sum(elem(ie)%mp(:,:)*elem(ie)%metdet(:,:))
301 : end do
302 1536 : call repro_sum(aratio, area, nelemd, nelemd, 1, commid=par%comm)
303 1536 : area(1) = 4.0_r8*pi/area(1) ! ratio correction
304 1536 : deallocate(aratio)
305 1536 : if (par%masterproc) then
306 2 : write(iulog,'(2a,f20.17)') subname, "re-initializing cube elements: area correction=", area(1)
307 2 : call shr_sys_flush(iulog)
308 : end if
309 :
310 12336 : do ie = 1, nelemd
311 10800 : call cube_init_atomic(elem(ie),gp%points,area(1))
312 12336 : call rotation_init_atomic(elem(ie),rot_type)
313 : end do
314 : end if
315 :
316 1536 : if(par%masterproc) then
317 2 : write(iulog,*) subname, 're-running mass_matrix'
318 2 : call shr_sys_flush(iulog)
319 : end if
320 1536 : call mass_matrix(par, elem)
321 :
322 : ! =================================================================
323 : ! Determine the global degree of freedome for each gridpoint
324 : ! =================================================================
325 1536 : if(par%masterproc) then
326 2 : write(iulog,*) subname, 'running global_dof'
327 2 : call shr_sys_flush(iulog)
328 : end if
329 1536 : call global_dof(par, elem)
330 :
331 : ! =================================================================
332 : ! Create Unique Indices
333 : ! =================================================================
334 :
335 12336 : do ie = 1, nelemd
336 12336 : call CreateUniqueIndex(elem(ie)%GlobalId,elem(ie)%gdofP,elem(ie)%idxP)
337 : end do
338 :
339 1536 : call SetElemOffset(par,elem, GlobalUniqueCols)
340 :
341 12336 : do ie = 1, nelemd
342 12336 : elem(ie)%idxV=>elem(ie)%idxP
343 : end do
344 :
345 : ! initialize flux terms to 0
346 12336 : do ie = 1, nelemd
347 43200000 : elem(ie)%derived%FM=0.0_r8
348 126630000 : elem(ie)%derived%FQ=0.0_r8
349 21103200 : elem(ie)%derived%FT=0.0_r8
350 21103200 : elem(ie)%derived%FDP=0.0_r8
351 21103200 : elem(ie)%derived%pecnd=0.0_r8
352 :
353 21103200 : elem(ie)%derived%Omega=0
354 63320400 : elem(ie)%state%dp3d=0
355 :
356 10800 : elem(ie)%derived%etadot_prescribed = nan
357 10800 : elem(ie)%derived%u_met = nan
358 10800 : elem(ie)%derived%v_met = nan
359 10800 : elem(ie)%derived%dudt_met = nan
360 10800 : elem(ie)%derived%dvdt_met = nan
361 10800 : elem(ie)%derived%T_met = nan
362 10800 : elem(ie)%derived%dTdt_met = nan
363 10800 : elem(ie)%derived%ps_met = nan
364 10800 : elem(ie)%derived%dpsdt_met = nan
365 10800 : elem(ie)%derived%nudge_factor = nan
366 :
367 17085600 : elem(ie)%derived%Utnd=0._r8
368 17085600 : elem(ie)%derived%Vtnd=0._r8
369 17087136 : elem(ie)%derived%Ttnd=0._r8
370 : end do
371 :
372 : ! ==========================================================
373 : ! This routines initalizes a Restart file. This involves:
374 : ! I) Setting up the MPI datastructures
375 : ! ==========================================================
376 1536 : deallocate(GridEdge)
377 8295936 : do j = 1, nelem
378 8295936 : call deallocate_gridvertex_nbrs(GridVertex(j))
379 : end do
380 1536 : deallocate(GridVertex)
381 :
382 12336 : do j = 1, MetaVertex(1)%nmembers
383 12336 : call deallocate_gridvertex_nbrs(MetaVertex(1)%members(j))
384 : end do
385 1536 : deallocate(MetaVertex)
386 :
387 : ! =====================================
388 : ! Set number of threads...
389 : ! =====================================
390 1536 : if(par%masterproc) then
391 2 : write(iulog,*) subname, "max_num_threads=",max_num_threads
392 2 : call shr_sys_flush(iulog)
393 : end if
394 :
395 1536 : nets = 1
396 1536 : nete = nelemd
397 1536 : call Prim_Advec_Init1(par, elem)
398 1536 : if (fv_nphys > 0) then
399 1536 : call fvm_init1(par,elem)
400 : end if
401 :
402 : ! =======================================================
403 : ! Allocate memory for subcell flux calculations.
404 : ! =======================================================
405 1536 : call allocate_subcell_integration_matrix_cslam(np, nc)
406 1536 : if (fv_nphys > 0) then
407 1536 : call allocate_subcell_integration_matrix_physgrid(np, fv_nphys)
408 : end if
409 :
410 1536 : call TimeLevel_init(tl)
411 :
412 1536 : if (fv_nphys > 0) then
413 1536 : if(par%masterproc) then
414 2 : write(iulog,*) subname, 'initialize basic fvm coordinate variables'
415 2 : call shr_sys_flush(iulog)
416 : end if
417 12336 : do ie = 1, nelemd
418 0 : call compute_basic_coordinate_vars(elem(ie), nc, irecons_tracer, &
419 0 : fvm(ie)%dalpha, fvm(ie)%dbeta, fvm(ie)%vtx_cart(:,:,1:nc,1:nc), &
420 0 : fvm(ie)%center_cart(1:nc,1:nc), fvm(ie)%area_sphere(1:nc,1:nc), &
421 10800 : fvm(ie)%spherecentroid(:,1:nc,1:nc))
422 0 : call compute_basic_coordinate_vars(elem(ie), fv_nphys, irecons_tracer,&
423 0 : fvm(ie)%dalpha_physgrid, fvm(ie)%dbeta_physgrid, &
424 21600 : fvm(ie)%vtx_cart_physgrid (:,:,1:fv_nphys,1:fv_nphys), &
425 21600 : fvm(ie)%center_cart_physgrid(1:fv_nphys,1:fv_nphys), &
426 21600 : fvm(ie)%area_sphere_physgrid(1:fv_nphys,1:fv_nphys), &
427 77136 : fvm(ie)%spherecentroid_physgrid(:,1:fv_nphys,1:fv_nphys))
428 : end do
429 : end if
430 :
431 1536 : if(par%masterproc) then
432 2 : write(iulog,*) subname, 'end of prim_init'
433 2 : call shr_sys_flush(iulog)
434 : end if
435 1536 : end subroutine prim_init1
436 : end module prim_init
|
