source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/domain.f90 @ 245

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

Creating temporary dynamico/lmdz/saturn branche

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