source: codes/icosagcm/trunk/src/domain.f90 @ 155

Last change on this file since 155 was 151, checked in by ymipsl, 11 years ago

Implementation of mixte parallelism MPI/OpenMP into src directory

YM
File size: 11.8 KB
Line 
1MODULE domain_mod
2  USE domain_param
3
4  TYPE t_domain
5    INTEGER :: face
6    INTEGER :: iim
7    INTEGER :: jjm
8    INTEGER :: ii_begin
9    INTEGER :: jj_begin
10    INTEGER :: ii_end
11    INTEGER :: jj_end
12    INTEGER :: ii_nb
13    INTEGER :: jj_nb
14    INTEGER :: ii_begin_glo
15    INTEGER :: jj_begin_glo
16    INTEGER :: ii_end_glo
17    INTEGER :: jj_end_glo
18    INTEGER,POINTER  :: assign_domain(:,:)
19    INTEGER,POINTER  :: assign_i(:,:)
20    INTEGER,POINTER  :: assign_j(:,:)
21    INTEGER,POINTER  :: edge_assign_domain(:,:,:)
22    INTEGER,POINTER  :: edge_assign_i(:,:,:)
23    INTEGER,POINTER  :: edge_assign_j(:,:,:)
24    INTEGER,POINTER  :: edge_assign_pos(:,:,:)
25    INTEGER,POINTER  :: edge_assign_sign(:,:,:)
26    REAL,POINTER     :: xyz(:,:,:)
27    REAL,POINTER     :: neighbour(:,:,:,:)
28    INTEGER,POINTER  :: delta(:,:)
29    REAL,POINTER     :: vertex(:,:,:,:)
30    LOGICAL,POINTER  :: own(:,:)
31    INTEGER,POINTER  :: ne(:,:,:)
32   
33    INTEGER :: t_right
34    INTEGER :: t_rup
35    INTEGER :: t_lup
36    INTEGER :: t_left
37    INTEGER :: t_ldown
38    INTEGER :: t_rdown
39
40    INTEGER :: u_right
41    INTEGER :: u_rup
42    INTEGER :: u_lup
43    INTEGER :: u_left
44    INTEGER :: u_ldown
45    INTEGER :: u_rdown
46
47    INTEGER :: z_rup
48    INTEGER :: z_up
49    INTEGER :: z_lup
50    INTEGER :: z_ldown
51    INTEGER :: z_down
52    INTEGER :: z_rdown
53   
54    INTEGER :: t_pos(6)
55    INTEGER :: u_pos(6)
56    INTEGER :: z_pos(6)
57     
58  END TYPE t_domain
59
60  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:)
61  INTEGER :: ndomain
62  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:)
63  INTEGER :: ndomain_glo
64
65  INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:)
66  INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:)
67  INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:)
68 
69CONTAINS
70
71  SUBROUTINE create_domain
72  USE grid_param
73  USE mpipara
74  IMPLICIT NONE
75  INTEGER :: ind,nf,ni,nj,i,j
76  INTEGER :: quotient, rest
77  TYPE(t_domain),POINTER :: d
78 
79    ndomain_glo=nsplit_i*nsplit_j*nb_face
80    ALLOCATE(domain_glo(ndomain_glo))
81    ALLOCATE(domglo_rank(ndomain_glo))
82    ALLOCATE(domglo_loc_ind(ndomain_glo))
83
84    ind=0
85    DO nf=1,nb_face
86      DO nj=1,nsplit_j
87        DO ni=1,nsplit_i
88          ind=ind+1
89          d=>domain_glo(ind)
90          quotient=(iim_glo/nsplit_i)
91          rest=MOD(iim_glo,nsplit_i)
92          IF (ni-1 < rest) THEN
93            d%ii_nb=quotient+1
94            d%ii_begin_glo=(quotient+1)*(ni-1)+1
95          ELSE
96            d%ii_nb=quotient
97            d%ii_begin_glo=(quotient+1)*rest+(ni-1-rest) * quotient+1
98          ENDIF
99          d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1
100 
101          IF (ni/=1) THEN
102            d%ii_nb=d%ii_nb+1
103            d%ii_begin_glo=d%ii_begin_glo-1
104          ENDIF
105         
106       
107          quotient=(jjm_glo/nsplit_j)
108          rest=MOD(jjm_glo,nsplit_j)
109          IF (nj-1 < rest) THEN
110            d%jj_nb=quotient+1
111            d%jj_begin_glo=(quotient+1)*(nj-1)+1
112          ELSE
113            d%jj_nb=quotient
114            d%jj_begin_glo=(quotient+1)*rest+(nj-1-rest) * quotient+1
115          ENDIF
116          d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1
117
118          IF (nj/=1) THEN
119            d%jj_nb=d%jj_nb+1
120            d%jj_begin_glo=d%jj_begin_glo-1
121          ENDIF
122
123
124          d%iim=d%ii_nb+2*halo
125          d%jjm=d%jj_nb+2*halo
126          d%ii_begin=halo+1
127          d%jj_begin=halo+1
128          d%ii_end=d%ii_begin+d%ii_nb-1
129          d%jj_end=d%jj_begin+d%jj_nb-1
130          d%face=nf       
131          ALLOCATE(d%assign_domain(d%iim,d%jjm))
132          ALLOCATE(d%assign_i(d%iim,d%jjm))
133          ALLOCATE(d%assign_j(d%iim,d%jjm))
134          ALLOCATE(d%edge_assign_domain(0:5,d%iim,d%jjm))
135          ALLOCATE(d%edge_assign_i(0:5,d%iim,d%jjm))
136          ALLOCATE(d%edge_assign_j(0:5,d%iim,d%jjm))
137          ALLOCATE(d%edge_assign_pos(0:5,d%iim,d%jjm))
138          ALLOCATE(d%edge_assign_sign(0:5,d%iim,d%jjm))
139          ALLOCATE(d%delta(d%iim,d%jjm))
140          ALLOCATE(d%xyz(3,d%iim,d%jjm))
141          ALLOCATE(d%vertex(3,0:5,d%iim,d%jjm))
142          ALLOCATE(d%neighbour(3,0:5,d%iim,d%jjm))
143          ALLOCATE(d%own(d%iim,d%jjm))
144          ALLOCATE(d%ne(0:5,d%iim,d%jjm))
145         
146          IF (is_mpi_root) PRINT *,"Domain ",ind," : size ",d%ii_nb,"x",d%jj_nb
147         
148        END DO
149      END DO
150    END DO
151   
152  END SUBROUTINE create_domain
153 
154  SUBROUTINE copy_domain(d1,d2)
155  IMPLICIT NONE
156  INTEGER :: face
157  TYPE(t_domain),TARGET,INTENT(IN) :: d1
158  TYPE(t_domain), INTENT(OUT) :: d2
159 
160    d2%iim = d1%iim
161    d2%jjm = d1%jjm
162    d2%ii_begin = d1%ii_begin
163    d2%jj_begin = d1%jj_begin
164    d2%ii_end = d1%ii_end
165    d2%jj_end = d1%jj_end
166    d2%ii_nb =  d1%ii_nb
167    d2%jj_nb = d1%jj_nb
168    d2%ii_begin_glo = d1%ii_begin_glo
169    d2%jj_begin_glo = d1%jj_begin_glo
170    d2%ii_end_glo = d1%ii_end_glo
171    d2%jj_end_glo = d1%jj_end_glo
172    d2%assign_domain => d1%assign_domain
173    d2%assign_i => d1%assign_i
174    d2%assign_j => d1%assign_j
175    d2%edge_assign_domain => d1%edge_assign_domain
176    d2%edge_assign_i => d1%edge_assign_i
177    d2%edge_assign_j => d1%edge_assign_j
178    d2%edge_assign_pos => d1%edge_assign_pos
179    d2%edge_assign_sign => d1%edge_assign_sign
180    d2%xyz => d1%xyz
181    d2%neighbour => d1%neighbour
182    d2%delta => d1%delta
183    d2%vertex => d1%vertex
184    d2%own => d1%own
185    d2%ne => d1%ne
186   
187    d2%t_right = d1%t_right
188    d2%t_rup = d1%t_rup
189    d2%t_lup = d1%t_lup
190    d2%t_left = d1%t_left
191    d2%t_ldown = d1%t_ldown
192    d2%t_rdown = d1%t_rdown
193
194    d2%u_right = d1%u_right
195    d2%u_rup = d1%u_rup
196    d2%u_lup = d1%u_lup
197    d2%u_left = d1%u_left
198    d2%u_ldown = d1%u_ldown
199    d2%u_rdown = d1%u_rdown
200
201    d2%z_rup = d1%z_rup
202    d2%z_up = d1%z_up
203    d2%z_lup = d1%z_lup
204    d2%z_ldown = d1%z_ldown
205    d2%z_down = d1%z_down
206    d2%z_rdown = d1%z_rdown
207   
208    d2%t_pos = d1%t_pos
209    d2%u_pos = d1%u_pos
210    d2%z_pos = d1%z_pos
211   
212  END SUBROUTINE copy_domain
213 
214 
215  SUBROUTINE assign_cell
216  USE metric
217  IMPLICIT NONE
218    INTEGER :: ind_d,ind,ind2,e
219    INTEGER :: nf,nf2
220    INTEGER :: i,j,k,ii,jj
221    TYPE(t_domain),POINTER :: d
222    INTEGER :: delta
223     
224   
225    DO ind_d=1,ndomain_glo
226      d=>domain_glo(ind_d)
227      nf=d%face
228      DO j=d%jj_begin,d%jj_end
229        DO i=d%ii_begin,d%ii_end
230          ii=d%ii_begin_glo-d%ii_begin+i
231          jj=d%jj_begin_glo-d%jj_begin+j
232          ind=vertex_glo(ii,jj,nf)%ind
233          delta=vertex_glo(ii,jj,nf)%delta
234          IF (cell_glo(ind)%assign_face==nf) THEN
235            cell_glo(ind)%assign_domain=ind_d
236            cell_glo(ind)%assign_i=i
237            cell_glo(ind)%assign_j=j
238          ENDIF
239          IF ( i+1 <= d%ii_end ) CALL assign_edge(ind_d,ind,i,j,delta,0)
240          IF ( j+1 <= d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,1)
241          IF ( i-1 >= d%ii_begin .AND. j+1<=d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,2)
242          IF ( i-1 >= d%ii_begin ) CALL assign_edge(ind_d,ind,i,j,delta,3)
243          IF ( j-1 >= d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,4)
244          IF ( i+1 <= d%ii_end .AND. j-1 >=d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,5)
245        ENDDO
246      ENDDO
247    ENDDO
248   
249   
250    DO ind_d=1,ndomain_glo
251      d=>domain_glo(ind_d)
252      nf=d%face
253      DO j=d%jj_begin-1,d%jj_end+1
254        DO i=d%ii_begin-1,d%ii_end+1
255          ii=d%ii_begin_glo-d%ii_begin+i
256          jj=d%jj_begin_glo-d%jj_begin+j
257          ind=vertex_glo(ii,jj,nf)%ind
258          d%assign_domain(i,j)=cell_glo(ind)%assign_domain
259          d%assign_i(i,j)=cell_glo(ind)%assign_i
260          d%assign_j(i,j)=cell_glo(ind)%assign_j
261          delta=vertex_glo(ii,jj,nf)%delta
262          d%delta(i,j)=vertex_glo(ii,jj,nf)%delta
263          DO k=0,5
264            ind2=vertex_glo(ii,jj,nf)%neighbour(k)
265            d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:)
266
267            d%ne(k,i,j)=1-2*MOD(k,2)
268
269            e=cell_glo(ind)%edge(MOD(k+delta+6,6))
270            d%edge_assign_domain(k,i,j)=edge_glo(e)%assign_domain
271            d%edge_assign_i(k,i,j)=edge_glo(e)%assign_i
272            d%edge_assign_j(k,i,j)=edge_glo(e)%assign_j
273            d%edge_assign_pos(k,i,j)=edge_glo(e)%assign_pos
274            nf2=domain_glo(edge_glo(e)%assign_domain)%face
275            d%edge_assign_sign(k,i,j)=1-2*MOD(12+tab_index(nf,nf2,0),2)
276            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) & 
277                /= MOD(edge_glo(e)%assign_pos+3,6)) THEN
278              d%edge_assign_sign(k,i,j)=-d%edge_assign_sign(k,i,j)
279            ENDIF
280             
281          ENDDO
282          d%xyz(:,i,j)=cell_glo(ind)%xyz(:)
283          IF (d%assign_domain(i,j)==ind_d) THEN
284           d%own(i,j)=.TRUE.
285          ELSE
286           d%own(i,j)=.FALSE.
287          ENDIF
288        ENDDO
289      ENDDO
290    ENDDO
291
292  CONTAINS
293
294    SUBROUTINE assign_edge(ind_d,ind,i,j,delta,k)
295      INTEGER :: ind_d,ind,i,j,delta,k
296      INTEGER :: e
297      e=cell_glo(ind)%edge(MOD(k+delta+6,6))
298      edge_glo(e)%assign_domain=ind_d
299      edge_glo(e)%assign_i=i
300      edge_glo(e)%assign_j=j
301      edge_glo(e)%assign_pos=k
302      edge_glo(e)%assign_delta=delta
303
304     END  SUBROUTINE assign_edge
305         
306  END SUBROUTINE assign_cell
307
308  SUBROUTINE compute_boundary
309  USE spherical_geom_mod
310  IMPLICIT NONE
311    INTEGER :: ind_d
312    INTEGER :: i,j,k
313    TYPE(t_domain),POINTER :: d 
314    REAL(rstd) :: ng1(3),ng2(3) 
315
316    DO ind_d=1,ndomain_glo
317      d=>domain_glo(ind_d)
318      DO j=d%jj_begin-1,d%jj_end+1
319        DO i=d%ii_begin-1,d%ii_end+1
320          DO k=0,5
321            ng1=d%neighbour(:,MOD(k,6),i,j)
322            ng2=d%neighbour(:,MOD(k+1,6),i,j)
323            IF (sqrt(sum((ng1-ng2)**2))<1e-16) ng2=d%neighbour(:,MOD(k+2,6),i,j)
324            CALL circumcenter(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j))
325          ENDDO
326        ENDDO
327      ENDDO
328    ENDDO       
329  END SUBROUTINE compute_boundary
330
331  SUBROUTINE set_neighbour_indice
332  USE metric
333  IMPLICIT NONE
334    INTEGER :: ind_d
335    TYPE(t_domain),POINTER :: d 
336   
337    DO ind_d=1,ndomain_glo
338      d=>domain_glo(ind_d)
339      d%t_right=1
340      d%t_left=-1
341      d%t_rup=d%iim
342      d%t_lup=d%iim-1
343      d%t_ldown=-d%iim
344      d%t_rdown=-d%iim+1
345     
346      d%u_right=0
347      d%u_lup=d%iim*d%jjm
348      d%u_ldown=2*d%iim*d%jjm
349     
350      d%u_rup=d%t_rup+d%u_ldown
351      d%u_left=d%t_left+d%u_right
352      d%u_rdown=d%t_rdown+d%u_lup
353     
354      d%z_up=0
355      d%z_down=d%iim*d%jjm
356      d%z_rup=d%t_rup+d%z_down
357      d%z_lup=d%t_lup+d%z_down
358      d%z_ldown=d%t_ldown+d%z_up
359      d%z_rdown=d%t_rdown+d%z_up
360     
361      d%t_pos(right)=d%t_right
362      d%t_pos(rup)=d%t_rup
363      d%t_pos(lup)=d%t_lup
364      d%t_pos(left)=d%t_left
365      d%t_pos(ldown)=d%t_ldown
366      d%t_pos(rdown)=d%t_rdown
367
368      d%u_pos(right)=d%u_right
369      d%u_pos(rup)=d%u_rup
370      d%u_pos(lup)=d%u_lup
371      d%u_pos(left)=d%u_left
372      d%u_pos(ldown)=d%u_ldown
373      d%u_pos(rdown)=d%u_rdown
374     
375      d%z_pos(vrup)=d%z_rup
376      d%z_pos(vup)=d%z_up
377      d%z_pos(vlup)=d%z_lup
378      d%z_pos(vldown)=d%z_ldown
379      d%z_pos(vdown)=d%z_down
380      d%z_pos(vrdown)=d%z_rdown
381     
382    ENDDO 
383 
384  END SUBROUTINE set_neighbour_indice
385     
386  SUBROUTINE assign_domain
387  USE mpipara
388  IMPLICIT NONE
389    INTEGER :: nb_domain(0:mpi_size-1)
390    INTEGER :: rank, ind,ind_glo
391   
392    DO rank=0,mpi_size-1
393      nb_domain(rank)=ndomain_glo/mpi_size
394      IF ( rank < MOD(ndomain_glo,mpi_size) ) nb_domain(rank)=nb_domain(rank)+1
395    ENDDO
396   
397    ndomain=nb_domain(mpi_rank)
398    ALLOCATE(domain(ndomain))
399    ALLOCATE(domloc_glo_ind(ndomain))
400   
401    rank=0
402    ind=0
403    DO ind_glo=1,ndomain_glo
404      ind=ind+1
405      domglo_rank(ind_glo)=rank
406      domglo_loc_ind(ind_glo)=ind
407      IF (rank==mpi_rank) THEN
408        CALL copy_domain(domain_glo(ind_glo),domain(ind))
409        domloc_glo_ind(ind)=ind_glo
410      ENDIF
411     
412      IF (ind==nb_domain(rank)) THEN
413        rank=rank+1
414        ind=0
415      ENDIF
416    ENDDO
417   
418  END SUBROUTINE  assign_domain     
419   
420         
421  SUBROUTINE compute_domain
422  IMPLICIT NONE
423    CALL init_domain_param
424    CALL create_domain
425    CALL assign_cell
426    CALL compute_boundary
427    CALL set_neighbour_indice
428    CALL assign_domain
429     
430  END SUBROUTINE compute_domain
431         
432END MODULE domain_mod 
433 
Note: See TracBrowser for help on using the repository browser.