1 | MODULE domain_mod |
---|
2 | USE domain_param |
---|
3 | IMPLICIT NONE |
---|
4 | SAVE |
---|
5 | |
---|
6 | TYPE t_cellset |
---|
7 | ! derived type containing arrays for lat,lon of center and bounds of cells of primal or dual mesh |
---|
8 | ! this information is needed to write CF-compliant NetCDF files |
---|
9 | ! and by XIOS |
---|
10 | INTEGER :: ncell |
---|
11 | INTEGER, ALLOCATABLE :: ij(:), & ! ij(n) = local index of cell n |
---|
12 | ind_glo(:), & ! ind_glo(n) = global index of cell n |
---|
13 | sgn(:) ! sign to apply before writing to disk (velocities) |
---|
14 | REAL, ALLOCATABLE :: lon(:), lat(:), bnds_lon(:,:), bnds_lat(:,:) |
---|
15 | END type t_cellset |
---|
16 | |
---|
17 | TYPE t_domain |
---|
18 | INTEGER :: face |
---|
19 | INTEGER :: iim |
---|
20 | INTEGER :: jjm |
---|
21 | INTEGER :: ii_begin |
---|
22 | INTEGER :: jj_begin |
---|
23 | INTEGER :: ii_end |
---|
24 | INTEGER :: jj_end |
---|
25 | INTEGER :: ii_nb |
---|
26 | INTEGER :: jj_nb |
---|
27 | INTEGER :: ii_begin_glo |
---|
28 | INTEGER :: jj_begin_glo |
---|
29 | INTEGER :: ii_end_glo |
---|
30 | INTEGER :: jj_end_glo |
---|
31 | INTEGER,POINTER :: assign_domain(:,:) |
---|
32 | INTEGER,POINTER :: assign_cell_glo(:,:) |
---|
33 | INTEGER,POINTER :: assign_i(:,:) |
---|
34 | INTEGER,POINTER :: assign_j(:,:) |
---|
35 | INTEGER,POINTER :: edge_assign_domain(:,:,:) |
---|
36 | INTEGER,POINTER :: edge_assign_i(:,:,:) |
---|
37 | INTEGER,POINTER :: edge_assign_j(:,:,:) |
---|
38 | INTEGER,POINTER :: edge_assign_pos(:,:,:) |
---|
39 | INTEGER,POINTER :: edge_assign_sign(:,:,:) |
---|
40 | INTEGER,POINTER :: vertex_assign_domain(:,:,:) |
---|
41 | INTEGER,POINTER :: vertex_assign_i(:,:,:) |
---|
42 | INTEGER,POINTER :: vertex_assign_j(:,:,:) |
---|
43 | INTEGER,POINTER :: vertex_assign_pos(:,:,:) |
---|
44 | REAL,POINTER :: xyz(:,:,:) |
---|
45 | REAL,POINTER :: neighbour(:,:,:,:) |
---|
46 | INTEGER,POINTER :: delta(:,:) |
---|
47 | REAL,POINTER :: vertex(:,:,:,:) |
---|
48 | LOGICAL,POINTER :: own(:,:) |
---|
49 | INTEGER,POINTER :: ne(:,:,:) |
---|
50 | |
---|
51 | INTEGER :: t_right |
---|
52 | INTEGER :: t_rup |
---|
53 | INTEGER :: t_lup |
---|
54 | INTEGER :: t_left |
---|
55 | INTEGER :: t_ldown |
---|
56 | INTEGER :: t_rdown |
---|
57 | |
---|
58 | INTEGER :: u_right |
---|
59 | INTEGER :: u_rup |
---|
60 | INTEGER :: u_lup |
---|
61 | INTEGER :: u_left |
---|
62 | INTEGER :: u_ldown |
---|
63 | INTEGER :: u_rdown |
---|
64 | |
---|
65 | INTEGER :: z_rup |
---|
66 | INTEGER :: z_up |
---|
67 | INTEGER :: z_lup |
---|
68 | INTEGER :: z_ldown |
---|
69 | INTEGER :: z_down |
---|
70 | INTEGER :: z_rdown |
---|
71 | |
---|
72 | INTEGER :: t_pos(6) |
---|
73 | INTEGER :: u_pos(6) |
---|
74 | INTEGER :: z_pos(6) |
---|
75 | |
---|
76 | TYPE(t_cellset), POINTER :: primal_own, dual_own, edge_own, & |
---|
77 | primal_all, dual_all, edge_all ! halo included, used for debug by writefield |
---|
78 | END TYPE t_domain |
---|
79 | |
---|
80 | TYPE t_mesh |
---|
81 | INTEGER :: ndomain, primal_num, edge_num, dual_num |
---|
82 | TYPE(t_domain), ALLOCATABLE :: domain(:) |
---|
83 | TYPE(t_cellset), ALLOCATABLE :: primal_own(:), dual_own(:), edge_own(:), & |
---|
84 | primal_all(:), dual_all(:), edge_all(:) |
---|
85 | END type t_mesh |
---|
86 | |
---|
87 | INTEGER :: ndomain |
---|
88 | INTEGER :: ndomain_glo |
---|
89 | |
---|
90 | LOGICAL :: swap_needed=.TRUE. ! .FALSE. if a thread always computes on the same domain |
---|
91 | !$OMP THREADPRIVATE(swap_needed) |
---|
92 | |
---|
93 | TYPE(t_mesh), TARGET :: mesh_loc, mesh_glo |
---|
94 | TYPE(t_domain), POINTER :: domain(:), domain_glo(:) |
---|
95 | |
---|
96 | INTEGER,ALLOCATABLE :: domglo_rank(:) ! size : ndomain_glo : mpi rank assigned to a domain |
---|
97 | INTEGER,ALLOCATABLE :: domglo_loc_ind(:) ! size : ndomain_glo : corresponding local indice for a global domain indice |
---|
98 | INTEGER,ALLOCATABLE :: domloc_glo_ind(:) ! size : ndomain : corresponding global indice for a local domain indice |
---|
99 | |
---|
100 | LOGICAL, ALLOCATABLE :: assigned_domain(:) ! size : ndomain : is an omp task is assigned to this domain |
---|
101 | !$OMP THREADPRIVATE(assigned_domain) |
---|
102 | |
---|
103 | CONTAINS |
---|
104 | |
---|
105 | SUBROUTINE allocate_mesh(ndom, mesh) |
---|
106 | INTEGER, INTENT(IN) :: ndom |
---|
107 | TYPE(t_mesh), TARGET, INTENT(OUT) :: mesh |
---|
108 | INTEGER :: ind |
---|
109 | mesh%ndomain = ndom |
---|
110 | ALLOCATE(mesh%domain(ndom)) |
---|
111 | ALLOCATE(mesh%primal_own(ndom)) |
---|
112 | ALLOCATE(mesh%primal_all(ndom)) |
---|
113 | ALLOCATE(mesh%dual_own(ndom)) |
---|
114 | ALLOCATE(mesh%dual_all(ndom)) |
---|
115 | ALLOCATE(mesh%edge_own(ndom)) |
---|
116 | ALLOCATE(mesh%edge_all(ndom)) |
---|
117 | DO ind=1, ndom |
---|
118 | mesh%domain(ind)%primal_own => mesh%primal_own(ind) |
---|
119 | mesh%domain(ind)%primal_all => mesh%primal_all(ind) |
---|
120 | mesh%domain(ind)%dual_own => mesh%dual_own(ind) |
---|
121 | mesh%domain(ind)%dual_all => mesh%dual_all(ind) |
---|
122 | mesh%domain(ind)%edge_own => mesh%edge_own(ind) |
---|
123 | mesh%domain(ind)%edge_all => mesh%edge_all(ind) |
---|
124 | END DO |
---|
125 | END SUBROUTINE allocate_mesh |
---|
126 | |
---|
127 | SUBROUTINE create_domain |
---|
128 | USE grid_param |
---|
129 | USE mpipara |
---|
130 | USE ioipsl |
---|
131 | INTEGER :: ind,nf,ni,nj,i,j |
---|
132 | INTEGER :: quotient, rest |
---|
133 | INTEGER :: halo_i,halo_j |
---|
134 | |
---|
135 | TYPE(t_domain),POINTER :: d |
---|
136 | |
---|
137 | ndomain_glo=nsplit_i*nsplit_j*nb_face |
---|
138 | CALL allocate_mesh(ndomain_glo, mesh_glo) |
---|
139 | domain_glo => mesh_glo%domain |
---|
140 | ALLOCATE(domglo_rank(ndomain_glo)) |
---|
141 | ALLOCATE(domglo_loc_ind(ndomain_glo)) |
---|
142 | |
---|
143 | halo_i=0 |
---|
144 | CALL getin("halo_i",halo_i) |
---|
145 | halo_j=1 |
---|
146 | CALL getin("halo_j",halo_j) |
---|
147 | |
---|
148 | ind=0 |
---|
149 | DO nf=1,nb_face |
---|
150 | DO nj=1,nsplit_j |
---|
151 | DO ni=1,nsplit_i |
---|
152 | ind=ind+1 |
---|
153 | d=>domain_glo(ind) |
---|
154 | quotient=(iim_glo/nsplit_i) |
---|
155 | rest=MOD(iim_glo,nsplit_i) |
---|
156 | IF (ni-1 < rest) THEN |
---|
157 | d%ii_nb=quotient+1 |
---|
158 | d%ii_begin_glo=(quotient+1)*(ni-1)+1 |
---|
159 | ELSE |
---|
160 | d%ii_nb=quotient |
---|
161 | d%ii_begin_glo=(quotient+1)*rest+(ni-1-rest) * quotient+1 |
---|
162 | ENDIF |
---|
163 | d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1 |
---|
164 | |
---|
165 | IF (ni/=1) THEN |
---|
166 | d%ii_nb=d%ii_nb+1 |
---|
167 | d%ii_begin_glo=d%ii_begin_glo-1 |
---|
168 | ENDIF |
---|
169 | |
---|
170 | |
---|
171 | quotient=(jjm_glo/nsplit_j) |
---|
172 | rest=MOD(jjm_glo,nsplit_j) |
---|
173 | IF (nj-1 < rest) THEN |
---|
174 | d%jj_nb=quotient+1 |
---|
175 | d%jj_begin_glo=(quotient+1)*(nj-1)+1 |
---|
176 | ELSE |
---|
177 | d%jj_nb=quotient |
---|
178 | d%jj_begin_glo=(quotient+1)*rest+(nj-1-rest) * quotient+1 |
---|
179 | ENDIF |
---|
180 | d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1 |
---|
181 | |
---|
182 | IF (nj/=1) THEN |
---|
183 | d%jj_nb=d%jj_nb+1 |
---|
184 | d%jj_begin_glo=d%jj_begin_glo-1 |
---|
185 | ENDIF |
---|
186 | |
---|
187 | |
---|
188 | d%iim=d%ii_nb+2*halo + halo_i*2 |
---|
189 | d%jjm=d%jj_nb+2*halo + halo_j*2 |
---|
190 | d%ii_begin=halo+1 + halo_i |
---|
191 | d%jj_begin=halo+1 + halo_j |
---|
192 | d%ii_end=d%ii_begin+d%ii_nb-1 |
---|
193 | d%jj_end=d%jj_begin+d%jj_nb-1 |
---|
194 | d%face=nf |
---|
195 | ALLOCATE(d%assign_domain(d%iim,d%jjm)) |
---|
196 | ALLOCATE(d%assign_cell_glo(d%iim,d%jjm)) |
---|
197 | ALLOCATE(d%assign_i(d%iim,d%jjm)) |
---|
198 | ALLOCATE(d%assign_j(d%iim,d%jjm)) |
---|
199 | ALLOCATE(d%edge_assign_domain(0:5,d%iim,d%jjm)) |
---|
200 | ALLOCATE(d%edge_assign_i(0:5,d%iim,d%jjm)) |
---|
201 | ALLOCATE(d%edge_assign_j(0:5,d%iim,d%jjm)) |
---|
202 | ALLOCATE(d%edge_assign_pos(0:5,d%iim,d%jjm)) |
---|
203 | ALLOCATE(d%edge_assign_sign(0:5,d%iim,d%jjm)) |
---|
204 | ALLOCATE(d%vertex_assign_domain(0:5,d%iim,d%jjm)) |
---|
205 | ALLOCATE(d%vertex_assign_i(0:5,d%iim,d%jjm)) |
---|
206 | ALLOCATE(d%vertex_assign_j(0:5,d%iim,d%jjm)) |
---|
207 | ALLOCATE(d%vertex_assign_pos(0:5,d%iim,d%jjm)) |
---|
208 | ALLOCATE(d%delta(d%iim,d%jjm)) |
---|
209 | ALLOCATE(d%xyz(3,d%iim,d%jjm)) |
---|
210 | ALLOCATE(d%vertex(3,0:5,d%iim,d%jjm)) |
---|
211 | ALLOCATE(d%neighbour(3,0:5,d%iim,d%jjm)) |
---|
212 | ALLOCATE(d%own(d%iim,d%jjm)) |
---|
213 | ALLOCATE(d%ne(0:5,d%iim,d%jjm)) |
---|
214 | |
---|
215 | IF (is_mpi_root) PRINT *,"Domain ",ind," : size ",d%ii_nb,"x",d%jj_nb |
---|
216 | |
---|
217 | END DO |
---|
218 | END DO |
---|
219 | END DO |
---|
220 | |
---|
221 | END SUBROUTINE create_domain |
---|
222 | |
---|
223 | SUBROUTINE copy_domain(d1,d2) |
---|
224 | INTEGER :: face |
---|
225 | TYPE(t_domain),TARGET,INTENT(IN) :: d1 |
---|
226 | TYPE(t_domain), INTENT(OUT) :: d2 |
---|
227 | |
---|
228 | d2%iim = d1%iim |
---|
229 | d2%jjm = d1%jjm |
---|
230 | d2%ii_begin = d1%ii_begin |
---|
231 | d2%jj_begin = d1%jj_begin |
---|
232 | d2%ii_end = d1%ii_end |
---|
233 | d2%jj_end = d1%jj_end |
---|
234 | d2%ii_nb = d1%ii_nb |
---|
235 | d2%jj_nb = d1%jj_nb |
---|
236 | d2%ii_begin_glo = d1%ii_begin_glo |
---|
237 | d2%jj_begin_glo = d1%jj_begin_glo |
---|
238 | d2%ii_end_glo = d1%ii_end_glo |
---|
239 | d2%jj_end_glo = d1%jj_end_glo |
---|
240 | d2%assign_domain => d1%assign_domain |
---|
241 | d2%assign_cell_glo => d1%assign_cell_glo |
---|
242 | d2%assign_i => d1%assign_i |
---|
243 | d2%assign_j => d1%assign_j |
---|
244 | d2%edge_assign_domain => d1%edge_assign_domain |
---|
245 | d2%edge_assign_i => d1%edge_assign_i |
---|
246 | d2%edge_assign_j => d1%edge_assign_j |
---|
247 | d2%edge_assign_pos => d1%edge_assign_pos |
---|
248 | d2%edge_assign_sign => d1%edge_assign_sign |
---|
249 | d2%vertex_assign_domain => d1%vertex_assign_domain |
---|
250 | d2%vertex_assign_i => d1%vertex_assign_i |
---|
251 | d2%vertex_assign_j => d1%vertex_assign_j |
---|
252 | d2%vertex_assign_pos => d1%vertex_assign_pos |
---|
253 | d2%xyz => d1%xyz |
---|
254 | d2%neighbour => d1%neighbour |
---|
255 | d2%delta => d1%delta |
---|
256 | d2%vertex => d1%vertex |
---|
257 | d2%own => d1%own |
---|
258 | d2%ne => d1%ne |
---|
259 | |
---|
260 | d2%t_right = d1%t_right |
---|
261 | d2%t_rup = d1%t_rup |
---|
262 | d2%t_lup = d1%t_lup |
---|
263 | d2%t_left = d1%t_left |
---|
264 | d2%t_ldown = d1%t_ldown |
---|
265 | d2%t_rdown = d1%t_rdown |
---|
266 | |
---|
267 | d2%u_right = d1%u_right |
---|
268 | d2%u_rup = d1%u_rup |
---|
269 | d2%u_lup = d1%u_lup |
---|
270 | d2%u_left = d1%u_left |
---|
271 | d2%u_ldown = d1%u_ldown |
---|
272 | d2%u_rdown = d1%u_rdown |
---|
273 | |
---|
274 | d2%z_rup = d1%z_rup |
---|
275 | d2%z_up = d1%z_up |
---|
276 | d2%z_lup = d1%z_lup |
---|
277 | d2%z_ldown = d1%z_ldown |
---|
278 | d2%z_down = d1%z_down |
---|
279 | d2%z_rdown = d1%z_rdown |
---|
280 | |
---|
281 | d2%t_pos = d1%t_pos |
---|
282 | d2%u_pos = d1%u_pos |
---|
283 | d2%z_pos = d1%z_pos |
---|
284 | |
---|
285 | END SUBROUTINE copy_domain |
---|
286 | |
---|
287 | |
---|
288 | SUBROUTINE assign_cell |
---|
289 | USE metric |
---|
290 | INTEGER :: ind_d,ind,ind2,e,v |
---|
291 | INTEGER :: nf,nf2 |
---|
292 | INTEGER :: i,j,k,ii,jj |
---|
293 | TYPE(t_domain),POINTER :: d |
---|
294 | INTEGER :: delta |
---|
295 | |
---|
296 | |
---|
297 | DO ind_d=1,ndomain_glo |
---|
298 | d=>domain_glo(ind_d) |
---|
299 | nf=d%face |
---|
300 | DO j=d%jj_begin,d%jj_end |
---|
301 | DO i=d%ii_begin,d%ii_end |
---|
302 | ii=d%ii_begin_glo-d%ii_begin+i |
---|
303 | jj=d%jj_begin_glo-d%jj_begin+j |
---|
304 | ind=vertex_glo(ii,jj,nf)%ind |
---|
305 | delta=vertex_glo(ii,jj,nf)%delta |
---|
306 | IF (cell_glo(ind)%assign_face==nf) THEN |
---|
307 | cell_glo(ind)%assign_domain=ind_d |
---|
308 | cell_glo(ind)%assign_i=i |
---|
309 | cell_glo(ind)%assign_j=j |
---|
310 | ENDIF |
---|
311 | IF ( i+1 <= d%ii_end ) CALL assign_edge(ind_d,ind,i,j,delta,0) |
---|
312 | IF ( j+1 <= d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,1) |
---|
313 | IF ( i-1 >= d%ii_begin .AND. j+1<=d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,2) |
---|
314 | IF ( i-1 >= d%ii_begin ) CALL assign_edge(ind_d,ind,i,j,delta,3) |
---|
315 | IF ( j-1 >= d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,4) |
---|
316 | IF ( i+1 <= d%ii_end .AND. j-1 >=d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,5) |
---|
317 | |
---|
318 | IF ( i+1 <= d%ii_end .AND. j+1 <= d%jj_end) CALL assign_vertex(ind_d,ind,i,j,delta,0) |
---|
319 | IF ( i-1 >= d%ii_begin .AND. j+1 <= d%jj_end) CALL assign_vertex(ind_d,ind,i,j,delta,1) |
---|
320 | IF ( i-1 >= d%ii_begin .AND. j+1 <= d%jj_end) CALL assign_vertex(ind_d,ind,i,j,delta,2) |
---|
321 | IF ( i-1 >= d%ii_begin .AND. j-1 >= d%jj_begin) CALL assign_vertex(ind_d,ind,i,j,delta,3) |
---|
322 | IF ( i+1 <= d%ii_end .AND. j-1 >= d%jj_begin) CALL assign_vertex(ind_d,ind,i,j,delta,4) |
---|
323 | IF ( i+1 <= d%ii_end .AND. j-1 >= d%jj_begin) CALL assign_vertex(ind_d,ind,i,j,delta,5) |
---|
324 | ENDDO |
---|
325 | ENDDO |
---|
326 | ENDDO |
---|
327 | |
---|
328 | |
---|
329 | DO ind_d=1,ndomain_glo |
---|
330 | d=>domain_glo(ind_d) |
---|
331 | nf=d%face |
---|
332 | DO j=d%jj_begin-1,d%jj_end+1 |
---|
333 | DO i=d%ii_begin-1,d%ii_end+1 |
---|
334 | ii=d%ii_begin_glo-d%ii_begin+i |
---|
335 | jj=d%jj_begin_glo-d%jj_begin+j |
---|
336 | ind=vertex_glo(ii,jj,nf)%ind |
---|
337 | d%assign_cell_glo(i,j) = ind |
---|
338 | d%assign_domain(i,j)=cell_glo(ind)%assign_domain |
---|
339 | d%assign_i(i,j)=cell_glo(ind)%assign_i |
---|
340 | d%assign_j(i,j)=cell_glo(ind)%assign_j |
---|
341 | delta=vertex_glo(ii,jj,nf)%delta |
---|
342 | d%delta(i,j)=vertex_glo(ii,jj,nf)%delta |
---|
343 | DO k=0,5 |
---|
344 | ind2=vertex_glo(ii,jj,nf)%neighbour(k) |
---|
345 | d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:) |
---|
346 | |
---|
347 | d%ne(k,i,j)=1-2*MOD(k,2) |
---|
348 | |
---|
349 | e=cell_glo(ind)%edge(MOD(k+delta+6,6)) |
---|
350 | d%edge_assign_domain(k,i,j)=edge_glo(e)%assign_domain |
---|
351 | d%edge_assign_i(k,i,j)=edge_glo(e)%assign_i |
---|
352 | d%edge_assign_j(k,i,j)=edge_glo(e)%assign_j |
---|
353 | d%edge_assign_pos(k,i,j)=edge_glo(e)%assign_pos |
---|
354 | nf2=domain_glo(edge_glo(e)%assign_domain)%face |
---|
355 | d%edge_assign_sign(k,i,j)=1-2*MOD(12+tab_index(nf,nf2,0),2) |
---|
356 | IF (MOD(6+k+tab_index(nf,nf2,0),6)/=edge_glo(e)%assign_pos .AND. MOD(6+k+tab_index(nf,nf2,0),6) & |
---|
357 | /= MOD(edge_glo(e)%assign_pos+3,6)) THEN |
---|
358 | d%edge_assign_sign(k,i,j)=-d%edge_assign_sign(k,i,j) |
---|
359 | ENDIF |
---|
360 | |
---|
361 | v=cell_glo(ind)%vertex(MOD(k+delta+6,6)) |
---|
362 | d%vertex_assign_domain(k,i,j)=vertices_glo(v)%assign_domain |
---|
363 | d%vertex_assign_i(k,i,j)=vertices_glo(v)%assign_i |
---|
364 | d%vertex_assign_j(k,i,j)=vertices_glo(v)%assign_j |
---|
365 | d%vertex_assign_pos(k,i,j)=vertices_glo(v)%assign_pos |
---|
366 | |
---|
367 | ENDDO |
---|
368 | d%xyz(:,i,j)=cell_glo(ind)%xyz(:) |
---|
369 | IF (d%assign_domain(i,j)==ind_d) THEN |
---|
370 | d%own(i,j)=.TRUE. |
---|
371 | ELSE |
---|
372 | d%own(i,j)=.FALSE. |
---|
373 | ENDIF |
---|
374 | ENDDO |
---|
375 | ENDDO |
---|
376 | ENDDO |
---|
377 | |
---|
378 | CONTAINS |
---|
379 | |
---|
380 | SUBROUTINE assign_edge(ind_d,ind,i,j,delta,k) |
---|
381 | INTEGER :: ind_d,ind,i,j,delta,k |
---|
382 | INTEGER :: e |
---|
383 | e=cell_glo(ind)%edge(MOD(k+delta+6,6)) |
---|
384 | edge_glo(e)%assign_domain=ind_d |
---|
385 | edge_glo(e)%assign_i=i |
---|
386 | edge_glo(e)%assign_j=j |
---|
387 | edge_glo(e)%assign_pos=k |
---|
388 | edge_glo(e)%assign_delta=delta |
---|
389 | |
---|
390 | END SUBROUTINE assign_edge |
---|
391 | |
---|
392 | SUBROUTINE assign_vertex(ind_d,ind,i,j,delta,k) |
---|
393 | INTEGER :: ind_d,ind,i,j,delta,k |
---|
394 | INTEGER :: e |
---|
395 | |
---|
396 | e=cell_glo(ind)%vertex(MOD(k+delta+6,6)) |
---|
397 | vertices_glo(e)%assign_domain=ind_d |
---|
398 | vertices_glo(e)%assign_i=i |
---|
399 | vertices_glo(e)%assign_j=j |
---|
400 | vertices_glo(e)%assign_pos=k |
---|
401 | vertices_glo(e)%assign_delta=delta |
---|
402 | |
---|
403 | END SUBROUTINE assign_vertex |
---|
404 | |
---|
405 | END SUBROUTINE assign_cell |
---|
406 | |
---|
407 | SUBROUTINE compute_boundary |
---|
408 | USE spherical_geom_mod |
---|
409 | INTEGER :: ind_d |
---|
410 | INTEGER :: i,j,k |
---|
411 | TYPE(t_domain),POINTER :: d |
---|
412 | REAL(rstd) :: ng1(3),ng2(3) |
---|
413 | |
---|
414 | DO ind_d=1,ndomain_glo |
---|
415 | d=>domain_glo(ind_d) |
---|
416 | DO j=d%jj_begin-1,d%jj_end+1 |
---|
417 | DO i=d%ii_begin-1,d%ii_end+1 |
---|
418 | DO k=0,5 |
---|
419 | ng1=d%neighbour(:,MOD(k,6),i,j) |
---|
420 | ng2=d%neighbour(:,MOD(k+1,6),i,j) |
---|
421 | IF (sqrt(sum((ng1-ng2)**2))<1e-16) ng2=d%neighbour(:,MOD(k+2,6),i,j) |
---|
422 | CALL circumcenter(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j)) |
---|
423 | ENDDO |
---|
424 | ENDDO |
---|
425 | ENDDO |
---|
426 | ENDDO |
---|
427 | END SUBROUTINE compute_boundary |
---|
428 | |
---|
429 | SUBROUTINE set_neighbour_indice |
---|
430 | USE metric |
---|
431 | INTEGER :: ind_d |
---|
432 | TYPE(t_domain),POINTER :: d |
---|
433 | |
---|
434 | DO ind_d=1,ndomain_glo |
---|
435 | d=>domain_glo(ind_d) |
---|
436 | d%t_right=1 |
---|
437 | d%t_left=-1 |
---|
438 | d%t_rup=d%iim |
---|
439 | d%t_lup=d%iim-1 |
---|
440 | d%t_ldown=-d%iim |
---|
441 | d%t_rdown=-d%iim+1 |
---|
442 | |
---|
443 | d%u_right=0 |
---|
444 | d%u_lup=d%iim*d%jjm |
---|
445 | d%u_ldown=2*d%iim*d%jjm |
---|
446 | |
---|
447 | d%u_rup=d%t_rup+d%u_ldown |
---|
448 | d%u_left=d%t_left+d%u_right |
---|
449 | d%u_rdown=d%t_rdown+d%u_lup |
---|
450 | |
---|
451 | d%z_up=0 |
---|
452 | d%z_down=d%iim*d%jjm |
---|
453 | d%z_rup=d%t_rup+d%z_down |
---|
454 | d%z_lup=d%t_lup+d%z_down |
---|
455 | d%z_ldown=d%t_ldown+d%z_up |
---|
456 | d%z_rdown=d%t_rdown+d%z_up |
---|
457 | |
---|
458 | d%t_pos(right)=d%t_right |
---|
459 | d%t_pos(rup)=d%t_rup |
---|
460 | d%t_pos(lup)=d%t_lup |
---|
461 | d%t_pos(left)=d%t_left |
---|
462 | d%t_pos(ldown)=d%t_ldown |
---|
463 | d%t_pos(rdown)=d%t_rdown |
---|
464 | |
---|
465 | d%u_pos(right)=d%u_right |
---|
466 | d%u_pos(rup)=d%u_rup |
---|
467 | d%u_pos(lup)=d%u_lup |
---|
468 | d%u_pos(left)=d%u_left |
---|
469 | d%u_pos(ldown)=d%u_ldown |
---|
470 | d%u_pos(rdown)=d%u_rdown |
---|
471 | |
---|
472 | d%z_pos(vrup)=d%z_rup |
---|
473 | d%z_pos(vup)=d%z_up |
---|
474 | d%z_pos(vlup)=d%z_lup |
---|
475 | d%z_pos(vldown)=d%z_ldown |
---|
476 | d%z_pos(vdown)=d%z_down |
---|
477 | d%z_pos(vrdown)=d%z_rdown |
---|
478 | |
---|
479 | ENDDO |
---|
480 | |
---|
481 | END SUBROUTINE set_neighbour_indice |
---|
482 | |
---|
483 | SUBROUTINE assign_domain |
---|
484 | USE mpipara |
---|
485 | USE grid_param |
---|
486 | INTEGER :: nb_domain(0:mpi_size-1) |
---|
487 | INTEGER :: rank, ind,ind_glo |
---|
488 | INTEGER :: block_j,jb,i,j,nd_glo,n,nf |
---|
489 | LOGICAL :: exit |
---|
490 | |
---|
491 | DO rank=0,mpi_size-1 |
---|
492 | nb_domain(rank)=ndomain_glo/mpi_size |
---|
493 | IF ( rank < MOD(ndomain_glo,mpi_size) ) nb_domain(rank)=nb_domain(rank)+1 |
---|
494 | ENDDO |
---|
495 | |
---|
496 | ndomain=nb_domain(mpi_rank) |
---|
497 | |
---|
498 | ! ALLOCATE(domain(ndomain)) |
---|
499 | CALL allocate_mesh(ndomain, mesh_loc) |
---|
500 | domain => mesh_loc%domain |
---|
501 | |
---|
502 | ALLOCATE(domloc_glo_ind(ndomain)) |
---|
503 | |
---|
504 | |
---|
505 | block_j=sqrt(nsplit_i*nsplit_j*nb_face*1./mpi_size) |
---|
506 | exit=.FALSE. |
---|
507 | jb=1 |
---|
508 | i=1 |
---|
509 | j=1 |
---|
510 | ind=1 |
---|
511 | nd_glo=0 |
---|
512 | rank=0 |
---|
513 | DO WHILE (.NOT. exit) |
---|
514 | |
---|
515 | IF (j==MIN(jb+block_j,nsplit_j*nb_face+1)) THEN |
---|
516 | j=jb |
---|
517 | i=i+1 |
---|
518 | ENDIF |
---|
519 | |
---|
520 | IF (i>nsplit_i) THEN |
---|
521 | i=1 |
---|
522 | jb=jb+block_j |
---|
523 | j=jb |
---|
524 | ENDIF |
---|
525 | |
---|
526 | IF (ind>nb_domain(rank)) THEN |
---|
527 | rank=rank+1 |
---|
528 | ind=1 |
---|
529 | ENDIF |
---|
530 | ind_glo=(j-1)*nsplit_i+i |
---|
531 | |
---|
532 | nd_glo=nd_glo+1 |
---|
533 | IF (nd_glo==ndomain_glo) THEN |
---|
534 | |
---|
535 | exit=.TRUE. |
---|
536 | IF (.NOT. (rank==mpi_size-1 .AND. ind==nb_domain(rank) )) THEN |
---|
537 | PRINT *, "Distribution problem in assign_domain" |
---|
538 | STOP |
---|
539 | ENDIF |
---|
540 | |
---|
541 | ENDIF |
---|
542 | |
---|
543 | domglo_rank(ind_glo)=rank |
---|
544 | domglo_loc_ind(ind_glo)=ind |
---|
545 | IF (rank==mpi_rank) THEN |
---|
546 | CALL copy_domain(domain_glo(ind_glo),domain(ind)) |
---|
547 | domloc_glo_ind(ind)=ind_glo |
---|
548 | ENDIF |
---|
549 | |
---|
550 | j=j+1 |
---|
551 | ind=ind+1 |
---|
552 | |
---|
553 | ENDDO |
---|
554 | |
---|
555 | IF (is_mpi_master) THEN |
---|
556 | |
---|
557 | ind_glo=0 |
---|
558 | PRINT *,'' |
---|
559 | PRINT*, ' MPI PROCESS DISTRIBUTION' |
---|
560 | PRINT *,'' |
---|
561 | |
---|
562 | WRITE(*,"(' ')", ADVANCE='NO') |
---|
563 | DO n=1,nsplit_i*7-1 |
---|
564 | WRITE(*,"('=')", ADVANCE='NO') |
---|
565 | ENDDO |
---|
566 | PRINT *,'' |
---|
567 | |
---|
568 | DO nf=1,nb_face |
---|
569 | DO j=1,nsplit_j |
---|
570 | IF (j>1) THEN |
---|
571 | WRITE(*,"(' ')", ADVANCE='NO') |
---|
572 | DO n=1,nsplit_i*7-1 |
---|
573 | WRITE(*,"('-')", ADVANCE='NO') |
---|
574 | ENDDO |
---|
575 | PRINT *,'' |
---|
576 | ENDIF |
---|
577 | |
---|
578 | WRITE(*,"('|')", ADVANCE='NO') |
---|
579 | DO i=1,nsplit_i |
---|
580 | WRITE(*,"(' ',' ',' |')",ADVANCE='NO') |
---|
581 | ENDDO |
---|
582 | PRINT *,'' |
---|
583 | |
---|
584 | WRITE(*,"('|')", ADVANCE='NO') |
---|
585 | DO i=1,nsplit_i |
---|
586 | ind_glo=ind_glo+1 |
---|
587 | WRITE(*,"(' ',i4.4 ,' |')", ADVANCE='NO') domglo_rank(ind_glo) |
---|
588 | END DO |
---|
589 | PRINT *,'' |
---|
590 | |
---|
591 | WRITE(*,"('|')", ADVANCE='NO') |
---|
592 | DO i=1,nsplit_i |
---|
593 | WRITE(*,"(' ',' ',' |')", ADVANCE='NO') |
---|
594 | ENDDO |
---|
595 | PRINT *,'' |
---|
596 | |
---|
597 | ENDDO |
---|
598 | |
---|
599 | WRITE(*,"(' ')", ADVANCE='NO') |
---|
600 | DO n=1,nsplit_i*7-1 |
---|
601 | WRITE(*,"('=')", ADVANCE='NO') |
---|
602 | ENDDO |
---|
603 | PRINT *,'' |
---|
604 | ENDDO |
---|
605 | ENDIF |
---|
606 | |
---|
607 | ! rank=0 |
---|
608 | ! ind=0 |
---|
609 | ! DO ind_glo=1,ndomain_glo |
---|
610 | ! ind=ind+1 |
---|
611 | ! domglo_rank(ind_glo)=rank |
---|
612 | ! domglo_loc_ind(ind_glo)=ind |
---|
613 | ! IF (rank==mpi_rank) THEN |
---|
614 | ! CALL copy_domain(domain_glo(ind_glo),domain(ind)) |
---|
615 | ! domloc_glo_ind(ind)=ind_glo |
---|
616 | ! ENDIF |
---|
617 | ! |
---|
618 | ! IF (ind==nb_domain(rank)) THEN |
---|
619 | ! rank=rank+1 |
---|
620 | ! ind=0 |
---|
621 | ! ENDIF |
---|
622 | ! ENDDO |
---|
623 | |
---|
624 | !$OMP PARALLEL |
---|
625 | CALL assign_domain_omp |
---|
626 | !$OMP END PARALLEL |
---|
627 | |
---|
628 | END SUBROUTINE assign_domain |
---|
629 | |
---|
630 | SUBROUTINE assign_domain_omp |
---|
631 | USE omp_para |
---|
632 | USE mpipara |
---|
633 | INTEGER :: nb_domain |
---|
634 | INTEGER :: rank, ind, i |
---|
635 | |
---|
636 | ALLOCATE(assigned_domain(ndomain)) |
---|
637 | |
---|
638 | ind=0 |
---|
639 | DO rank=0,omp_domain_size-1 |
---|
640 | nb_domain=ndomain/omp_domain_size |
---|
641 | IF ( rank < MOD(ndomain,omp_domain_size) ) nb_domain=nb_domain+1 |
---|
642 | |
---|
643 | DO i=1,nb_domain |
---|
644 | ind=ind+1 |
---|
645 | IF (rank==omp_domain_rank) THEN |
---|
646 | assigned_domain(ind)=.TRUE. |
---|
647 | WRITE(*,'(A,I6.6,A,I4.4,A,I4.4,A,I8.8)') "Rank ",mpi_rank," task ",rank," local domain ",ind, & |
---|
648 | " global domain ",domloc_glo_ind(ind) |
---|
649 | ELSE |
---|
650 | assigned_domain(ind)=.FALSE. |
---|
651 | ENDIF |
---|
652 | ENDDO |
---|
653 | |
---|
654 | ENDDO |
---|
655 | |
---|
656 | END SUBROUTINE assign_domain_omp |
---|
657 | |
---|
658 | SUBROUTINE init_domain_unst(dom) |
---|
659 | USE grid_param, ONLY : primal_num |
---|
660 | TYPE(t_domain) :: dom |
---|
661 | dom%ii_begin=1 |
---|
662 | dom%ii_end=primal_num |
---|
663 | dom%iim=primal_num |
---|
664 | dom%jj_begin=1 |
---|
665 | dom%jj_end=1 |
---|
666 | dom%jjm=1 |
---|
667 | ALLOCATE(dom%assign_domain(primal_num,1)) |
---|
668 | dom%assign_domain(:,:)=1 ! all cells are assigned to MPI rank 0 |
---|
669 | END SUBROUTINE init_domain_unst |
---|
670 | |
---|
671 | SUBROUTINE compute_domain |
---|
672 | USE grid_param, ONLY : grid_type, grid_unst, grid_ico |
---|
673 | INTEGER :: ind |
---|
674 | SELECT CASE(grid_type) |
---|
675 | CASE(grid_unst) |
---|
676 | ! FIXME temporary, sequential hack |
---|
677 | ndomain_glo=1 |
---|
678 | CALL allocate_mesh(ndomain_glo, mesh_glo) |
---|
679 | domain_glo => mesh_glo%domain |
---|
680 | ALLOCATE(domglo_rank(ndomain_glo)) |
---|
681 | ALLOCATE(domglo_loc_ind(ndomain_glo)) |
---|
682 | domglo_rank(:)=0 |
---|
683 | domglo_loc_ind(:)=1 |
---|
684 | |
---|
685 | ndomain=1 |
---|
686 | CALL allocate_mesh(ndomain, mesh_loc) |
---|
687 | domain => mesh_loc%domain |
---|
688 | ALLOCATE(domloc_glo_ind(ndomain)) |
---|
689 | ALLOCATE(assigned_domain(ndomain)) |
---|
690 | domloc_glo_ind(:)=1 |
---|
691 | assigned_domain(:)=.TRUE. |
---|
692 | |
---|
693 | CALL init_domain_unst(domain(1)) |
---|
694 | CALL init_domain_unst(domain_glo(1)) |
---|
695 | |
---|
696 | CASE DEFAULT |
---|
697 | CALL init_domain_param |
---|
698 | CALL create_domain |
---|
699 | CALL assign_cell |
---|
700 | CALL compute_boundary |
---|
701 | CALL set_neighbour_indice |
---|
702 | CALL assign_domain |
---|
703 | END SELECT |
---|
704 | |
---|
705 | ! ALLOCATE(primal_cells(ndomain)) |
---|
706 | ! ALLOCATE(dual_cells(ndomain)) |
---|
707 | ! ALLOCATE(primal_cells_single(ndomain)) |
---|
708 | ! ALLOCATE(dual_cells_single(ndomain)) |
---|
709 | ! DO ind=1, ndomain |
---|
710 | ! domain(ind)%primal_cells => primal_cells(ind) |
---|
711 | ! domain(ind)%dual_cells => dual_cells(ind) |
---|
712 | ! domain(ind)%primal_cells_single => primal_cells_single(ind) |
---|
713 | ! domain(ind)%dual_cells_single => dual_cells_single(ind) |
---|
714 | ! END DO |
---|
715 | |
---|
716 | END SUBROUTINE compute_domain |
---|
717 | |
---|
718 | END MODULE domain_mod |
---|