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

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

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

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