source: XIOS3/trunk/src/test/test_ugrid.f90 @ 2433

Last change on this file since 2433 was 2409, checked in by jderouillat, 21 months ago

Backport commits [2181,2200-2202,2250,2252] related to UGRID in XIOS3
File size: 9.8 KB
Line 
1PROGRAM test_ugrid
2
3  USE xios
4  USE mod_wait
5  IMPLICIT NONE
6  INCLUDE "mpif.h"
7  INTEGER :: comm_rank
8  INTEGER :: comm_size
9  INTEGER :: ierr
10
11  CHARACTER(len=*),PARAMETER :: id="client"
12  INTEGER :: comm
13  TYPE(xios_duration) :: dtime
14  CHARACTER(len=20) :: dtime_str
15  TYPE(xios_date) :: date
16  CHARACTER(len=20) :: date_str
17  CHARACTER(len=15) :: calendar_type
18  TYPE(xios_context) :: ctx_hdl
19  INTEGER,PARAMETER :: ni_glo=100
20  INTEGER,PARAMETER :: nj_glo=100
21  INTEGER,PARAMETER :: llm=5
22  DOUBLE PRECISION  :: lval(llm)=1, scalar = 5
23  TYPE(xios_field) :: field_hdl
24  TYPE(xios_fieldgroup) :: fieldgroup_hdl
25  TYPE(xios_file) :: file_hdl
26  LOGICAL :: ok
27
28  DOUBLE PRECISION,DIMENSION(ni_glo,nj_glo) :: lon_glo,lat_glo
29  DOUBLE PRECISION :: field_A_glo(ni_glo,nj_glo,llm)
30  DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:),field_A(:,:,:), lonvalue(:,:), axisValue(:), field_domain(:,:) ;
31  INTEGER :: ni,ibegin,iend,nj,jbegin,jend
32  INTEGER :: i,j,l,ts,n
33
34  REAL,PARAMETER    :: Pi = 4.*atan(1.)
35
36  INTEGER,PARAMETER :: ncells=6
37  REAL :: cellx(ncells),celly(ncells),cellz(ncells)
38  REAL :: cellx_bounds(4,ncells),celly_bounds(4,ncells),cellz_bounds(4,ncells)
39  REAL(8) :: cell_lon(ncells), cell_lat(ncells)
40  REAL(8) :: cell_bounds_lon(4,ncells), cell_bounds_lat(4,ncells)
41  REAL :: field_cell(ncells) = (/1,2,3,4,5,6/)
42
43  INTEGER,PARAMETER :: nedges=12
44  REAL :: edgex(nedges),edgey(nedges),edgez(nedges)
45  REAL :: edgex_bounds(2,nedges),edgey_bounds(2,nedges),edgez_bounds(2,nedges)
46  REAL(8) :: edge_lon(nedges), edge_lat(nedges)
47  REAL(8) :: edge_bounds_lon(2,nedges), edge_bounds_lat(2,nedges)
48  REAL :: field_edge(nedges) = (/1,2,3,4,5,6,7,8,9,10,11,12/)
49
50  INTEGER,PARAMETER :: nnodes=8
51  REAL :: nodex(nnodes),nodey(nnodes),nodez(nnodes)
52  REAL(8) :: node_lon(nnodes), node_lat(nnodes)
53  REAL :: field_node(nnodes) = (/1,2,3,4,5,6,7,8/)
54 
55   
56  cellx(1) = 1
57  celly(1) = 0 
58  cellz(1) = 0 
59  cellx_bounds(1,1)=  1 ; cellx_bounds(2,1)=  1 ; cellx_bounds(3,1)=  1  ; cellx_bounds(4,1)=  1
60  celly_bounds(1,1)= -1 ; celly_bounds(2,1)=  1 ; celly_bounds(3,1)=  1  ; celly_bounds(4,1)= -1
61  cellz_bounds(1,1)= -1 ; cellz_bounds(2,1)= -1 ; cellz_bounds(3,1)=  1  ; cellz_bounds(4,1)=  1
62 
63  cellx(2) = 0
64  celly(2) = 1 
65  cellz(2) = 0 
66  cellx_bounds(1,2)= -1 ; cellx_bounds(2,2)=  1 ; cellx_bounds(3,2)=  1  ; cellx_bounds(4,2)= -1
67  celly_bounds(1,2)=  1 ; celly_bounds(2,2)=  1 ; celly_bounds(3,2)=  1  ; celly_bounds(4,2)=  1
68  cellz_bounds(1,2)= -1 ; cellz_bounds(2,2)= -1 ; cellz_bounds(3,2)=  1  ; cellz_bounds(4,2)=  1
69   
70  cellx(3) = -1
71  celly(3) = 0 
72  cellz(3) = 0 
73  cellx_bounds(1,3)= -1 ; cellx_bounds(2,3)= -1 ; cellx_bounds(3,3)= -1  ; cellx_bounds(4,3)= -1
74  celly_bounds(1,3)= -1 ; celly_bounds(2,3)=  1 ; celly_bounds(3,3)=  1  ; celly_bounds(4,3)= -1
75  cellz_bounds(1,3)= -1 ; cellz_bounds(2,3)= -1 ; cellz_bounds(3,3)=  1  ; cellz_bounds(4,3)=  1
76 
77  cellx(4) = 0
78  celly(4) = -1 
79  cellz(4) = 0 
80  cellx_bounds(1,4)= -1 ; cellx_bounds(2,4)=  1 ; cellx_bounds(3,4)=  1  ; cellx_bounds(4,4)= -1
81  celly_bounds(1,4)= -1 ; celly_bounds(2,4)= -1 ; celly_bounds(3,4)= -1  ; celly_bounds(4,4)= -1
82  cellz_bounds(1,4)= -1 ; cellz_bounds(2,4)= -1 ; cellz_bounds(3,4)=  1  ; cellz_bounds(4,4)=  1
83   
84  cellx(5) = 0
85  celly(5) = 0 
86  cellz(5) = 1 
87  cellx_bounds(1,5)= -1 ; cellx_bounds(2,5)=  1 ; cellx_bounds(3,5)=  1  ; cellx_bounds(4,5)= -1
88  celly_bounds(1,5)= -1 ; celly_bounds(2,5)= -1 ; celly_bounds(3,5)=  1  ; celly_bounds(4,5)=  1
89  cellz_bounds(1,5)=  1 ; cellz_bounds(2,5)=  1 ; cellz_bounds(3,5)=  1  ; cellz_bounds(4,5)=  1
90 
91  cellx(6) = 0
92  celly(6) = 0 
93  cellz(6) = -1 
94  cellx_bounds(1,6)= -1 ; cellx_bounds(2,6)=  1 ; cellx_bounds(3,6)=  1  ; cellx_bounds(4,6)= -1
95  celly_bounds(1,6)= -1 ; celly_bounds(2,6)= -1 ; celly_bounds(3,6)=  1  ; celly_bounds(4,6)=  1
96  cellz_bounds(1,6)= -1 ; cellz_bounds(2,6)= -1 ; cellz_bounds(3,6)= -1  ; cellz_bounds(4,6)= -1
97
98  CALL cube2lonlat_1(cellx,celly,cellz,cell_lon,cell_lat)
99  CALL cube2lonlat_2(cellx_bounds,celly_bounds,cellz_bounds,cell_bounds_lon,cell_bounds_lat)
100
101 
102  edgex(1) =  1
103  edgey(1) =  0
104  edgez(1) =  1
105  edgex_bounds(1,1) =  1 ; edgex_bounds(2,1) =  1
106  edgey_bounds(1,1) = -1 ; edgey_bounds(2,1) =  1
107  edgez_bounds(1,1) =  1 ; edgez_bounds(2,1) =  1
108 
109  edgex(2) = -1
110  edgey(2) =  0
111  edgez(2) =  1
112  edgex_bounds(1,2) = -1 ; edgex_bounds(2,2) = -1
113  edgey_bounds(1,2) = -1 ; edgey_bounds(2,2) =  1
114  edgez_bounds(1,2) =  1 ; edgez_bounds(2,2) =  1
115 
116  edgex(3) =  0
117  edgey(3) =  1
118  edgez(3) =  1
119  edgex_bounds(1,3) = -1 ; edgex_bounds(2,3) =  1
120  edgey_bounds(1,3) =  1 ; edgey_bounds(2,3) =  1
121  edgez_bounds(1,3) =  1 ; edgez_bounds(2,3) =  1
122 
123  edgex(4) =  0
124  edgey(4) = -1
125  edgez(4) =  1
126  edgex_bounds(1,4) = -1 ; edgex_bounds(2,4) =  1
127  edgey_bounds(1,4) = -1 ; edgey_bounds(2,4) = -1
128  edgez_bounds(1,4) =  1 ; edgez_bounds(2,4) =  1
129
130  edgex(5) =  1
131  edgey(5) =  0
132  edgez(5) = -1
133  edgex_bounds(1,5) =  1 ; edgex_bounds(2,5) =  1
134  edgey_bounds(1,5) = -1 ; edgey_bounds(2,5) =  1
135  edgez_bounds(1,5) = -1 ; edgez_bounds(2,5) = -1
136 
137  edgex(6) = -1
138  edgey(6) =  0
139  edgez(6) =  -1
140  edgex_bounds(1,6) = -1 ; edgex_bounds(2,6) = -1
141  edgey_bounds(1,6) = -1 ; edgey_bounds(2,6) =  1
142  edgez_bounds(1,6) = -1 ; edgez_bounds(2,6) = -1
143 
144  edgex(7) =  0
145  edgey(7) =  1
146  edgez(7) =  -1
147  edgex_bounds(1,7) = -1 ; edgex_bounds(2,7) =  1
148  edgey_bounds(1,7) =  1 ; edgey_bounds(2,7) =  1
149  edgez_bounds(1,7) = -1 ; edgez_bounds(2,7) = -1
150 
151  edgex(8) =  0
152  edgey(8) = -1
153  edgez(8) =  -1
154  edgex_bounds(1,8) = -1 ; edgex_bounds(2,8) =  1
155  edgey_bounds(1,8) = -1 ; edgey_bounds(2,8) = -1
156  edgez_bounds(1,8) = -1 ; edgez_bounds(2,8) = -1
157 
158  edgex(9) =  1
159  edgey(9) =  1
160  edgez(9) =  0
161  edgex_bounds(1,9) =  1 ; edgex_bounds(2,9) =  1
162  edgey_bounds(1,9) =  1 ; edgey_bounds(2,9) =  1
163  edgez_bounds(1,9) = -1 ; edgez_bounds(2,9) =  1
164
165  edgex(10) = -1
166  edgey(10) =  1
167  edgez(10) =  0
168  edgex_bounds(1,10) = -1 ; edgex_bounds(2,10) = -1
169  edgey_bounds(1,10) =  1 ; edgey_bounds(2,10) =  1
170  edgez_bounds(1,10) = -1 ; edgez_bounds(2,10) =  1
171
172  edgex(11) = -1
173  edgey(11) = -1
174  edgez(11) =  0
175  edgex_bounds(1,11) = -1 ; edgex_bounds(2,11) = -1
176  edgey_bounds(1,11) = -1 ; edgey_bounds(2,11) = -1
177  edgez_bounds(1,11) = -1 ; edgez_bounds(2,11) =  1
178
179  edgex(12) =  1
180  edgey(12) = -1
181  edgez(12) =  0
182  edgex_bounds(1,12) =  1 ; edgex_bounds(2,12) =  1
183  edgey_bounds(1,12) = -1 ; edgey_bounds(2,12) = -1
184  edgez_bounds(1,12) = -1 ; edgez_bounds(2,12) =  1
185
186  CALL cube2lonlat_1(edgex, edgey, edgez, edge_lon, edge_lat)
187  CALL cube2lonlat_2(edgex_bounds, edgey_bounds, edgez_bounds, edge_bounds_lon, edge_bounds_lat)
188
189  nodex(1) = -1
190  nodey(1) = -1
191  nodez(1) = -1
192
193  nodex(2) = 1
194  nodey(2) = -1
195  nodez(2) = -1
196
197  nodex(3) = 1
198  nodey(3) = 1
199  nodez(3) = -1
200
201  nodex(4) = -1
202  nodey(4) = 1
203  nodez(4) = -1
204
205  nodex(5) = -1
206  nodey(5) = -1
207  nodez(5) =  1
208
209  nodex(6) =  1
210  nodey(6) = -1
211  nodez(6) =  1
212
213  nodex(7) =  1
214  nodey(7) =  1
215  nodez(7) =  1
216
217  nodex(8) = -1
218  nodey(8) =  1
219  nodez(8) =  1
220
221  CALL cube2lonlat_1(nodex,nodey,nodez,node_lon,node_lat)
222   
223!!! MPI Initialization
224
225  CALL MPI_INIT(ierr)
226
227  CALL init_wait
228
229!!! XIOS Initialization (get the local communicator)
230
231  CALL xios_initialize(id,return_comm=comm)
232  CALL xios_context_initialize(id, comm)
233 
234  CALL MPI_COMM_RANK(comm,comm_rank,ierr)
235  CALL MPI_COMM_SIZE(comm,comm_size,ierr)
236
237  CALL xios_set_domain_attr("cells",ni_glo=ncells, ibegin=0, ni=ncells, nvertex=4 , type='unstructured')
238  CALL xios_set_domain_attr("cells",lonvalue_1d=cell_lon, latvalue_1d=cell_lat, bounds_lon_1d=cell_bounds_lon, bounds_lat_1d= cell_bounds_lat)
239 
240  CALL xios_set_domain_attr("edges",ni_glo=nedges, ibegin=0, ni=nedges, nvertex=2 , type='unstructured')
241  CALL xios_set_domain_attr("edges",lonvalue_1d=edge_lon, latvalue_1d=edge_lat, bounds_lon_1d=edge_bounds_lon, bounds_lat_1d= edge_bounds_lat)
242 
243  CALL xios_set_domain_attr("nodes",ni_glo=nnodes, ibegin=0, ni=nnodes, nvertex=1 , type='unstructured')
244  CALL xios_set_domain_attr("nodes",lonvalue_1d=node_lon, latvalue_1d=node_lat)
245
246  dtime%second = 3600
247  CALL xios_set_timestep(dtime)
248  CALL xios_close_context_definition
249 
250  DO ts=1,1
251    CALL xios_update_calendar(ts)
252    CALL xios_send_field("field_node",field_node) 
253    CALL xios_send_field("field_edge",field_edge) 
254    CALL xios_send_field("field_cell",field_cell) 
255    CALL wait_us(5000) ;
256  ENDDO
257
258  CALL xios_context_finalize()
259  CALL xios_finalize()
260
261  CALL MPI_FINALIZE(ierr)
262
263CONTAINS
264 
265  SUBROUTINE cube2sphere(xc,yc,zc,xs,ys,zs)
266  IMPLICIT NONE
267  REAL,INTENT(IN)  :: xc, yc, zc
268  REAL,INTENT(OUT) :: xs, ys, zs
269 
270    xs=xc/sqrt(1-yc*yc/2-zc*zc/2+yc*yc*zc*zc/3)
271    ys=yc/sqrt(1-xc*xc/2-zc*zc/2+xc*xc*zc*zc/3)
272    zs=zc/sqrt(1-xc*xc/2-yc*yc/2+xc*xc*yc*yc/3)
273  END SUBROUTINE cube2sphere
274 
275  SUBROUTINE xyz2lonlat(x,y,z,lon,lat)
276  IMPLICIT NONE
277  REAL,INTENT(IN)  :: x,y,z
278  REAL(8),INTENT(OUT) :: lon,lat
279  REAL             :: norm
280 
281    norm=sqrt(x*x+y*y+z*z)
282 
283    lat=asin(z/norm)*180./Pi
284    lon=atan2(y/norm,x/norm)*180/Pi
285  END SUBROUTINE xyz2lonlat
286 
287  SUBROUTINE cube2lonlat_0(xc,yc,zc,lon,lat)
288  IMPLICIT NONE
289    REAL,INTENT(IN)  :: xc,yc,zc
290    REAL(8),INTENT(OUT) :: lon,lat
291    REAL             :: xs,ys,zs
292    CALL cube2sphere(xc,yc,zc,xs,ys,zs)
293    CALL xyz2lonlat(xs,ys,zs,lon,lat)
294  END SUBROUTINE  cube2lonlat_0
295
296  SUBROUTINE cube2lonlat_1(xc,yc,zc,lon,lat)
297  IMPLICIT NONE
298  REAL,INTENT(IN),DIMENSION(:)  :: xc,yc,zc
299  REAL(8),INTENT(OUT),DIMENSION(:) :: lon,lat
300  INTEGER :: i
301    DO i=1,SIZE(xc,1)
302      CALL cube2lonlat_0(xc(i),yc(i),zc(i),lon(i),lat(i))
303    ENDDO
304  END SUBROUTINE cube2lonlat_1
305 
306  SUBROUTINE cube2lonlat_2(xc,yc,zc,lon,lat)
307  IMPLICIT NONE
308  REAL,INTENT(IN),DIMENSION(:,:)  :: xc,yc,zc
309  REAL(8),INTENT(OUT),DIMENSION(:,:) :: lon,lat
310  INTEGER :: i,j
311 
312    DO j=1,SIZE(xc,2)
313      DO i=1,SIZE(xc,1)
314        CALL cube2lonlat_0(xc(i,j),yc(i,j),zc(i,j),lon(i,j),lat(i,j))
315      ENDDO
316    ENDDO
317
318  END SUBROUTINE cube2lonlat_2
319
320END PROGRAM test_ugrid
321
322
323
324
325
Note: See TracBrowser for help on using the repository browser.