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

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