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

Last change on this file since 266 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

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