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

Last change on this file since 50 was 26, checked in by ymipsl, 12 years ago

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