Context Navigation

source: XIOS/dev/dev_trunk_graph/src/test/generic_testcase.f90 @ 2019

Last change on this file since 2019 was 2019, checked in by yushan, 3 years ago
Graph intermedia commit to a tmp branch
File size: 87.2 KB

Line
1	PROGRAM generic_testcase
2	USE xios
3	USE mod_wait
4	IMPLICIT NONE
5	INCLUDE "mpif.h"
6	INTEGER,PARAMETER :: unit=10
7	INTEGER,PARAMETER :: len_str=255
8	INTEGER :: ierr
9
10	INTEGER :: nb_proc_atm
11	INTEGER :: nb_proc_oce
12	INTEGER :: nb_proc_surf
13	CHARACTER(LEN=len_str) :: proc_distribution
14	CHARACTER(LEN=len_str) :: start_date
15	CHARACTER(LEN=len_str) :: duration
16	DOUBLE PRECISION :: sypd
17	NAMELIST /params_run/ nb_proc_atm, nb_proc_surf, nb_proc_oce, &
18	duration, sypd, start_date, proc_distribution
19
20
21	TYPE tmodel_params
22	CHARACTER(len_str) :: timestep=""
23	CHARACTER(len_str) :: domain=""
24	DOUBLE PRECISION :: domain_proc_frac
25	INTEGER :: domain_proc_n
26	CHARACTER(len_str) :: axis=""
27	DOUBLE PRECISION :: axis_proc_frac
28	INTEGER :: axis_proc_n
29	INTEGER :: ensemble_proc_n
30	INTEGER :: ni
31	INTEGER :: nj
32	INTEGER :: nlev
33	CHARACTER(len_str) :: init_field2D=""
34	DOUBLE PRECISION :: pressure_factor
35	LOGICAL :: domain_mask
36	CHARACTER(len_str) :: domain_mask_type=""
37	LOGICAL :: scalar_mask
38	CHARACTER(len_str) :: scalar_mask_type=""
39	LOGICAL :: axis_mask
40	LOGICAL :: mask3d
41	INTEGER :: field_sub_freq
42	INTEGER :: field_sub_offset
43	END TYPE tmodel_params
44
45
46
47	CHARACTER(len_str) :: timestep_atm
48	CHARACTER(len_str) :: domain_atm
49	CHARACTER(len_str) :: axis_atm
50	INTEGER :: ni_atm
51	INTEGER :: nj_atm
52	INTEGER :: nlev_atm
53	CHARACTER(len_str) :: init_field2D_atm
54
55
56	LOGICAL :: i_am_atm, i_am_oce, i_am_surf, i_am_server
57	INTEGER :: rank, size_loc
58
59
60	OPEN(unit=unit, file='param.def',status='old',iostat=ierr)
61	nb_proc_atm=1
62	nb_proc_oce=0
63	nb_proc_surf=0
64	duration="1d"
65	sypd=10000
66	start_date="2018-01-01"
67	proc_distribution="contiguous"
68	READ (unit, nml=params_run)
69	CLOSE(unit)
70	!!! MPI Initialization
71
72	CALL MPI_INIT(ierr)
73	CALL init_wait
74	CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
75	CALL MPI_COMM_SIZE(MPI_COMM_WORLD,size_loc,ierr)
76	if (size_loc < nb_proc_atm + nb_proc_oce + nb_proc_surf) then
77	PRINT *,"inconsistency : size=",size_loc," nb_proc_atm=",nb_proc_atm,&
78	" nb_proc_oce=",nb_proc_oce," nb_proc_surf=",nb_proc_surf
79	CALL MPI_ABORT()
80	endif
81
82	i_am_atm=.FALSE. ; i_am_oce=.FALSE. ; i_am_surf=.FALSE. ; i_am_server=.FALSE.
83
84	IF (proc_distribution=="contiguous") THEN
85
86	! First ranks are dedicated to Atm, next to Ocean, next to Surf, last to Server
87
88	IF (rank < nb_proc_atm) THEN
89	i_am_atm=.TRUE.
90	ELSE IF (rank < nb_proc_atm+nb_proc_oce) THEN
91	i_am_oce=.TRUE.
92	ELSE IF (rank < nb_proc_atm+nb_proc_oce+nb_proc_surf) THEN
93	i_am_surf=.TRUE.
94	ELSE
95	i_am_server=.TRUE.
96	ENDIF ;
97
98	ENDIF
99
100
101	IF (i_am_server) THEN
102	CALL xios_init_server()
103
104	ELSE
105
106	IF (i_am_atm) CALL model("atm")
107	IF (i_am_oce) CALL model("oce")
108	IF (i_am_surf) CALL model("surf")
109
110	CALL xios_finalize()
111	ENDIF
112
113	CALL MPI_FINALIZE(ierr)
114
115	if(i_am_server) then
116	print*, "server ",rank, " finished"
117	else
118	print*, "client ",rank, " finished"
119	endif
120
121	CONTAINS
122
123	SUBROUTINE model(model_id)
124	IMPLICIT NONE
125	CHARACTER(LEN=*) :: model_id
126	TYPE(tmodel_params) :: params, other_params
127	TYPE(xios_context) :: ctx_hdl
128	INTEGER :: comm
129	TYPE(xios_date) :: cdate, edate
130	TYPE(xios_duration) :: dtime
131	INTEGER :: rank
132	INTEGER :: size_loc
133	DOUBLE PRECISION :: timestep_in_seconds, simulated_seconds_per_seconds, elapsed_per_timestep
134	INTEGER :: ts
135
136	DOUBLE PRECISION, POINTER :: lon(:)
137	DOUBLE PRECISION, POINTER :: lat(:)
138	DOUBLE PRECISION, POINTER :: axis_value(:)
139	LOGICAL, POINTER :: domain_mask(:)
140	INTEGER, POINTER :: domain_index(:)
141	LOGICAL, POINTER :: axis_mask(:)
142	LOGICAL :: scalar_mask
143	INTEGER, POINTER :: axis_index(:)
144	DOUBLE PRECISION, POINTER :: ensemble_value(:)
145
146	DOUBLE PRECISION, POINTER :: X_lon(:)
147	DOUBLE PRECISION, POINTER :: X_lat(:)
148	LOGICAL, POINTER :: X_mask(:)
149	INTEGER, POINTER :: X_index(:)
150
151	DOUBLE PRECISION, POINTER :: Y_lon(:)
152	DOUBLE PRECISION, POINTER :: Y_lat(:)
153	LOGICAL, POINTER :: Y_mask(:)
154	INTEGER, POINTER :: Y_index(:)
155
156	DOUBLE PRECISION, POINTER :: field2D_init(:)
157	DOUBLE PRECISION, POINTER :: fieldX_init(:)
158	DOUBLE PRECISION, POINTER :: fieldY_init(:)
159	DOUBLE PRECISION, POINTER :: fieldXY_init(:,:)
160
161	DOUBLE PRECISION, POINTER :: field2D(:), other_field2D(:)
162	DOUBLE PRECISION :: field0D, other_field0D
163	DOUBLE PRECISION, POINTER :: field_X(:), other_field_X(:)
164	DOUBLE PRECISION, POINTER :: field_Y(:), other_field_Y(:)
165	DOUBLE PRECISION, POINTER :: field_Z(:), other_field_Z(:)
166	DOUBLE PRECISION, POINTER :: field_XY(:,:), other_field_XY(:,:)
167	DOUBLE PRECISION, POINTER :: field_XYZ(:,:,:), other_field_XYZ(:,:,:)
168	DOUBLE PRECISION, POINTER :: field_XZ(:,:), other_field_XZ(:,:)
169	DOUBLE PRECISION, POINTER :: field_YZ(:,:), other_field_YZ(:,:)
170
171	DOUBLE PRECISION, POINTER :: field2D_sub(:), other_field2D_sub(:)
172	DOUBLE PRECISION, POINTER :: field3D(:,:), other_field3D(:,:)
173	DOUBLE PRECISION, POINTER :: field3D_sub(:,:), other_field3D_sub(:,:)
174	DOUBLE PRECISION, POINTER :: field3D_recv(:,:), other_field3D_recv(:,:)
175	DOUBLE PRECISION, POINTER :: pressure(:,:), other_pressure(:,:)
176
177	DOUBLE PRECISION, POINTER :: field2D_W(:,:), other_field2D_W(:,:)
178	DOUBLE PRECISION, POINTER :: field0D_W(:), other_field0D_W(:)
179	DOUBLE PRECISION, POINTER :: field_XW(:,:), other_field_XW(:,:)
180	DOUBLE PRECISION, POINTER :: field_YW(:,:), other_field_YW(:,:)
181	DOUBLE PRECISION, POINTER :: field_ZW(:,:), other_field_ZW(:,:)
182	DOUBLE PRECISION, POINTER :: field_XYW(:,:,:), other_field_XYW(:,:,:)
183	DOUBLE PRECISION, POINTER :: field_XYZW(:,:,:,:), other_field_XYZW(:,:,:,:)
184	DOUBLE PRECISION, POINTER :: field_XZW(:,:,:), other_field_XZW(:,:,:)
185	DOUBLE PRECISION, POINTER :: field_YZW(:,:,:), other_field_YZW(:,:,:)
186
187	DOUBLE PRECISION, POINTER :: field2D_sub_W(:,:), other_field2D_sub_W(:,:)
188	DOUBLE PRECISION, POINTER :: field3D_W(:,:,:), other_field3D_W(:,:,:)
189	DOUBLE PRECISION, POINTER :: field3D_sub_W(:,:,:), other_field3D_sub_W(:,:,:)
190	DOUBLE PRECISION, POINTER :: field3D_recv_W(:,:,:), other_field3D_recv_W(:,:,:)
191	DOUBLE PRECISION, POINTER :: pressure_W(:,:,:), other_pressure_W(:,:,:)
192
193
194
195
196	TYPE(xios_duration) :: freq_op_recv, freq_offset_recv
197	INTEGER :: ni,nj,nk
198	INTEGER :: i,j,k,xy,x,y,z,w
199	DOUBLE PRECISION :: scale,dist
200	LOGICAL :: ok
201	INTEGER :: ierr
202
203	LOGICAL :: ok_field0D, ok_field2D, ok_field3D, ok_pressure, ok_field2D_sub, ok_field3D_sub,ok_field3D_recv, ok_field3D_send
204	LOGICAL :: ok_field_X, ok_field_Y, ok_field_XY, ok_field_Z, ok_field_XYZ, ok_field_XZ, ok_field_YZ
205	LOGICAL :: ok_field0D_W, ok_field2D_W, ok_field3D_W, ok_pressure_W, ok_field2D_sub_W, ok_field3D_sub_W,ok_field3D_recv_W, ok_field3D_send_W
206	LOGICAL :: ok_field_XW, ok_field_YW, ok_field_XYW, ok_field_ZW, ok_field_XYZW, ok_field_XZW, ok_field_YZW
207
208	LOGICAL :: ok_other_field0D, ok_other_field2D, ok_other_field3D, ok_other_pressure, ok_other_field2D_sub, ok_other_field3D_sub,ok_other_field3D_recv, ok_other_field3D_send
209	LOGICAL :: ok_other_field_X, ok_other_field_Y, ok_other_field_XY, ok_other_field_Z, ok_other_field_XYZ, ok_other_field_XZ, ok_other_field_YZ
210	LOGICAL :: ok_other_field0D_W, ok_other_field2D_W, ok_other_field3D_W, ok_other_pressure_W, ok_other_field2D_sub_W, ok_other_field3D_sub_W,ok_other_field3D_recv_W, ok_other_field3D_send_W
211	LOGICAL :: ok_other_field_XW, ok_other_field_YW, ok_other_field_XYW, ok_other_field_ZW, ok_other_field_XYZW, ok_other_field_XZW, ok_other_field_YZW
212
213	!! XIOS Initialization (get the local communicator)
214	CALL xios_initialize(trim(model_id),return_comm=comm)
215	CALL MPI_COMM_RANK(comm,rank,ierr)
216	CALL MPI_COMM_SIZE(comm,size_loc,ierr)
217	CALL xios_context_initialize(trim(model_id),comm)
218	CALL xios_get_handle(trim(model_id),ctx_hdl)
219	CALL xios_set_current_context(ctx_hdl)
220
221	CALL init_model_params('',params)
222	CALL init_model_params('other_',other_params)
223	!!! Definition du dÃ©coupage domain/axis
224
225
226	!!! Definition de la start date et du timestep
227	CALL xios_set_start_date(xios_date_convert_from_string(start_date//" 00:00:00"))
228	dtime=xios_duration_convert_from_string(TRIM(params%timestep))
229	CALL xios_set_timestep(timestep=dtime)
230
231	!!! Calcul de temps elaps par seconde pour respecter le SYPD (hyp : pas de delai d'I/O)
232	CALL xios_get_start_date(cdate)
233	edate=cdate+xios_duration_convert_from_string(TRIM(duration))
234	timestep_in_seconds=xios_date_convert_to_seconds(cdate+dtime) - xios_date_convert_to_seconds(cdate)
235	simulated_seconds_per_seconds=sypd * 365
236	elapsed_per_timestep=timestep_in_seconds/simulated_seconds_per_seconds ! in seconds
237
238
239	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
240	!!! standard initialisation of domain, axis, grid, field !
241	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242
243	CALL init_domain("domain", comm, params, ni, nj, lon, lat, domain_mask, domain_index, &
244	X_lon, X_lat, X_mask, X_index, &
245	Y_lon, Y_lat, Y_mask, Y_index)
246
247	CALL init_axis("axis", comm, params, axis_value, axis_mask, axis_index)
248	CALL init_scalar("scalar", comm, params, scalar_mask)
249	CALL init_ensemble("ensemble", comm, params, ensemble_value)
250
251	CALL set_mask3d("grid3D",params, ni, nj, lon, lat , axis_value)
252	!!! Fin de la definition du contexte
253
254	ok_field3D_recv=xios_is_valid_field("field3D_recv").AND.xios_is_valid_field("field3D_resend") ;
255	IF (ok_field3D_recv) THEN
256	CALL xios_is_defined_field_attr("field3D_recv",freq_op=ok)
257	IF (ok) THEN
258	CALL xios_get_field_attr("field3D_recv",freq_op=freq_op_recv)
259	ELSE
260	freq_op_recv%timestep=1
261	CALL xios_set_field_attr("field3D_recv",freq_op=freq_op_recv)
262	ENDIF
263	CALL xios_is_defined_field_attr("field3D_recv",freq_offset=ok)
264	IF (ok) THEN
265	CALL xios_get_field_attr("field3D_recv",freq_offset=freq_offset_recv)
266	ELSE
267	freq_offset_recv%timestep=0
268	CALL xios_set_field_attr("field3D_recv",freq_op=freq_op_recv)
269	ENDIF
270	CALL xios_set_field_attr("field3D_resend",freq_op=freq_op_recv,freq_offset=freq_offset_recv)
271	ENDIF
272
273	ok_field3D_recv_W=xios_is_valid_field("field3D_recv_W").AND.xios_is_valid_field("field3D_resend_W") ;
274	IF (ok_field3D_recv_W) THEN
275	CALL xios_is_defined_field_attr("field3D_recv_W",freq_op=ok)
276	IF (ok) THEN
277	CALL xios_get_field_attr("field3D_recv_W",freq_op=freq_op_recv)
278	ELSE
279	freq_op_recv%timestep=1
280	CALL xios_set_field_attr("field3D_recv_W",freq_op=freq_op_recv)
281	ENDIF
282	CALL xios_is_defined_field_attr("field3D_recv_W",freq_offset=ok)
283	IF (ok) THEN
284	CALL xios_get_field_attr("field3D_recv_W",freq_offset=freq_offset_recv)
285	ELSE
286	freq_offset_recv%timestep=0
287	CALL xios_set_field_attr("field3D_recv_W",freq_op=freq_op_recv)
288	ENDIF
289	CALL xios_set_field_attr("field3D_resend_W",freq_op=freq_op_recv,freq_offset=freq_offset_recv)
290	ENDIF
291
292
293	CALL init_field2D(comm, params, lon, lat, domain_mask, field2D_init, &
294	x_lon, x_lat, x_mask, fieldX_init, y_lon, y_lat, y_mask, fieldY_init, fieldXY_init)
295
296	xy=size(domain_index)
297	x=size(X_index)
298	y=size(Y_index)
299	z=size(axis_index)
300	w=size(ensemble_value)
301
302	ALLOCATE(field2D(0:xy-1))
303	ALLOCATE(field3D(0:xy-1,0:z-1))
304	ALLOCATE(pressure(0:xy-1,0:z-1))
305	ALLOCATE(field3D_recv(0:xy-1,0:z-1))
306	ALLOCATE(field_Z(0:z-1))
307	ALLOCATE(field_X(0:x-1))
308	ALLOCATE(field_Y(0:y-1))
309	ALLOCATE(field_XY(0:x-1,0:y-1))
310	ALLOCATE(field_XYZ(0:x-1,0:y-1,0:z-1))
311	ALLOCATE(field_XZ(0:x-1,0:z-1))
312	ALLOCATE(field_YZ(0:y-1,0:z-1))
313
314	ALLOCATE(field2D_W(0:xy-1,0:w-1))
315	ALLOCATE(field3D_W(0:xy-1,0:z-1,0:w-1))
316	ALLOCATE(pressure_W(0:xy-1,0:z-1,0:w-1))
317	ALLOCATE(field3D_recv_W(0:xy-1,0:z-1,0:w-1))
318	ALLOCATE(field_ZW(0:z-1,0:w-1))
319	ALLOCATE(field_XW(0:x-1,0:w-1))
320	ALLOCATE(field_YW(0:y-1,0:w-1))
321	ALLOCATE(field_XYW(0:x-1,0:y-1,0:w-1))
322	ALLOCATE(field_XYZW(0:x-1,0:y-1,0:z-1,0:w-1))
323	ALLOCATE(field_XZW(0:x-1,0:z-1,0:w-1))
324	ALLOCATE(field_YZW(0:y-1,0:z-1,0:w-1))
325	ALLOCATE(field0D_W(0:w-1))
326
327
328
329	DO i=0,xy-1
330	IF (domain_index(i)/=-1) THEN
331	field2D(i)=field2D_init(domain_index(i))
332	ENDIF
333	ENDDO
334
335	DO i=0,x-1
336	IF (X_index(i)/=-1) THEN
337	field_X(i)=fieldX_init(X_index(i))
338	ENDIF
339	ENDDO
340
341	DO i=0,y-1
342	IF (Y_index(i)/=-1) THEN
343	field_Y(i)=fieldY_init(Y_index(i))
344	ENDIF
345	ENDDO
346
347	DO j=0,y-1
348	DO i=0,x-1
349	IF (Y_index(j)/=-1 .AND. X_index(i)/=-1) THEN
350	field_XY(i,j)=fieldXY_init(X_index(i),Y_index(j))
351	ENDIF
352	ENDDO
353	ENDDO
354
355	DO i=0,z-1
356	IF (axis_index(i)/=-1) THEN
357	k=axis_index(i)
358	IF (axis_mask(k)) THEN
359	field_Z(i)=axis_value(axis_index(k))
360	field3D(:,i)=field2D(:)*axis_value(axis_index(k))
361	field_XYZ(:,:,i)=field_XY(:,:)*axis_value(axis_index(k))
362	field_XZ(:,i)=field_X(:)*axis_value(axis_index(k))
363	field_YZ(:,i)=field_Y(:)*axis_value(axis_index(k))
364	ENDIF
365	ENDIF
366	ENDDO
367
368	pressure=1e20
369	DO j=0,z-1
370	pressure(:,j)=axis_value(j) ;
371	DO i=0,xy-1
372	IF (domain_index(i)/=-1) THEN
373	k=domain_index(i)
374	IF (domain_mask(k)) THEN
375	dist=sqrt((lat(k)/90.)2+(lon(k)/180.)2) ;
376	pressure(i,j)=pressure(i,j)(1+params%pressure_factorexp(-(4dist)*2))
377	ENDIF
378	ENDIF
379	ENDDO
380	ENDDO
381
382	field0D=1
383
384
385	field2D_W(:,0) = field2D(:)(1+0.1ensemble_value(0))
386	field3D_W(:,:,0)= field3D(:,:)(1+0.1ensemble_value(0))
387	pressure_W(:,:,0) = pressure(:,:)(1+0.1ensemble_value(0))
388	field_XW(:,0) = field_X(:)(1+0.1ensemble_value(0))
389	field_YW(:,0) = field_Y(:)(1+0.1ensemble_value(0))
390	field_XYW(:,:,0) = field_XY(:,:)(1+0.1ensemble_value(0))
391	field_ZW(:,0) = field_Z(:)(1+0.1ensemble_value(0))
392	field_XYZW(:,:,:,0) = field_XYZ(:,:,:)(1+0.1ensemble_value(0))
393	field_XZW(:,:,0) = field_XZ(:,:)(1+0.1ensemble_value(0))
394	field_YZW(:,:,0) = field_YZ(:,:)(1+0.1ensemble_value(0))
395	field0D_W(0) = field0D(1+0.1ensemble_value(0))
396
397	ok_field2D=xios_is_valid_field("field2D") ;
398	ok_field3D=xios_is_valid_field("field3D") ;
399	ok_pressure=xios_is_valid_field("pressure") ;
400	ok_field2D_sub=xios_is_valid_field("field2D_sub") ;
401	ok_field3D_sub=xios_is_valid_field("field3D_sub") ;
402	ok_field_X=xios_is_valid_field("field_X") ;
403	ok_field_Y=xios_is_valid_field("field_Y") ;
404	ok_field_XY=xios_is_valid_field("field_XY") ;
405	ok_field_Z=xios_is_valid_field("field_Z") ;
406	ok_field_XYZ=xios_is_valid_field("field_XYZ") ;
407	ok_field_XZ=xios_is_valid_field("field_XZ") ;
408	ok_field_YZ=xios_is_valid_field("field_YZ") ;
409	ok_field0D=xios_is_valid_field("field0D") ;
410
411	ok_field2D_W=xios_is_valid_field("field2D_W") ;
412	ok_field3D_W=xios_is_valid_field("field3D_W") ;
413	ok_pressure_W=xios_is_valid_field("pressure_W") ;
414	ok_field2D_sub_W=xios_is_valid_field("field2D_sub_W") ;
415	ok_field3D_sub_W=xios_is_valid_field("field3D_sub_W") ;
416	ok_field_XW=xios_is_valid_field("field_XW") ;
417	ok_field_YW=xios_is_valid_field("field_YW") ;
418	ok_field_XYW=xios_is_valid_field("field_XYW") ;
419	ok_field_ZW=xios_is_valid_field("field_ZW") ;
420	ok_field_XYZW=xios_is_valid_field("field_XYZW") ;
421	ok_field_XZW=xios_is_valid_field("field_XZW") ;
422	ok_field_YZW=xios_is_valid_field("field_YZW") ;
423	ok_field0D_W=xios_is_valid_field("field0D_W") ;
424
425
426
427	IF (ASSOCIATED(lon)) DEALLOCATE(lon)
428	IF (ASSOCIATED(lat)) DEALLOCATE(lat)
429	IF (ASSOCIATED(axis_value)) DEALLOCATE(axis_value)
430	IF (ASSOCIATED(domain_mask)) DEALLOCATE(domain_mask)
431	IF (ASSOCIATED(domain_index)) DEALLOCATE(domain_index)
432	IF (ASSOCIATED(axis_mask)) DEALLOCATE(axis_mask)
433	IF (ASSOCIATED(axis_index)) DEALLOCATE(axis_index)
434	IF (ASSOCIATED(ensemble_value)) DEALLOCATE(ensemble_value)
435	IF (ASSOCIATED(X_lon)) DEALLOCATE(X_lon)
436	IF (ASSOCIATED(X_lat)) DEALLOCATE(X_lat)
437	IF (ASSOCIATED(X_mask)) DEALLOCATE(X_mask)
438	IF (ASSOCIATED(X_index)) DEALLOCATE(X_index)
439	IF (ASSOCIATED(Y_lon)) DEALLOCATE(Y_lon)
440	IF (ASSOCIATED(Y_lat)) DEALLOCATE(Y_lat)
441	IF (ASSOCIATED(Y_mask)) DEALLOCATE(Y_mask)
442	IF (ASSOCIATED(Y_index)) DEALLOCATE(Y_index)
443	IF (ASSOCIATED(field2D_init)) DEALLOCATE(field2D_init)
444	IF (ASSOCIATED(fieldX_init)) DEALLOCATE(fieldX_init)
445	IF (ASSOCIATED(fieldY_init)) DEALLOCATE(fieldY_init)
446	IF (ASSOCIATED(fieldXY_init)) DEALLOCATE(fieldXY_init)
447
448
449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450	! duplicated section for other_ !
451	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452
453	CALL init_domain("other_domain", comm, other_params, ni, nj, lon, lat, domain_mask, domain_index, &
454	X_lon, X_lat, X_mask, X_index, &
455	Y_lon, Y_lat, Y_mask, Y_index)
456
457	CALL init_axis("other_axis", comm, other_params, axis_value, axis_mask, axis_index)
458	CALL init_scalar("other_scalar", comm, params, scalar_mask)
459	CALL init_ensemble("other_ensemble", comm, other_params, ensemble_value)
460
461	CALL set_mask3d("other_grid3D",other_params, ni, nj, lon, lat , axis_value)
462	!!! Fin de la definition du contexte
463
464	ok_other_field3D_recv=xios_is_valid_field("other_field3D_recv").AND.xios_is_valid_field("other_field3D_resend") ;
465	IF (ok_other_field3D_recv) THEN
466	CALL xios_is_defined_field_attr("other_field3D_recv",freq_op=ok)
467	IF (ok) THEN
468	CALL xios_get_field_attr("other_field3D_recv",freq_op=freq_op_recv)
469	ELSE
470	freq_op_recv%timestep=1
471	CALL xios_set_field_attr("other_field3D_recv",freq_op=freq_op_recv)
472	ENDIF
473	CALL xios_is_defined_field_attr("other_field3D_recv",freq_offset=ok)
474	IF (ok) THEN
475	CALL xios_get_field_attr("other_field3D_recv",freq_offset=freq_offset_recv)
476	ELSE
477	freq_offset_recv%timestep=0
478	CALL xios_set_field_attr("other_field3D_recv",freq_op=freq_op_recv)
479	ENDIF
480	CALL xios_set_field_attr("other_field3D_resend",freq_op=freq_op_recv,freq_offset=freq_offset_recv)
481	ENDIF
482
483	ok_other_field3D_recv_W=xios_is_valid_field("other_field3D_recv_W").AND.xios_is_valid_field("other_field3D_resend_W") ;
484	IF (ok_other_field3D_recv_W) THEN
485	CALL xios_is_defined_field_attr("other_field3D_recv_W",freq_op=ok)
486	IF (ok) THEN
487	CALL xios_get_field_attr("other_field3D_recv_W",freq_op=freq_op_recv)
488	ELSE
489	freq_op_recv%timestep=1
490	CALL xios_set_field_attr("other_field3D_recv_W",freq_op=freq_op_recv)
491	ENDIF
492	CALL xios_is_defined_field_attr("other_field3D_recv_W",freq_offset=ok)
493	IF (ok) THEN
494	CALL xios_get_field_attr("other_field3D_recv_W",freq_offset=freq_offset_recv)
495	ELSE
496	freq_offset_recv%timestep=0
497	CALL xios_set_field_attr("other_field3D_recv_W",freq_op=freq_op_recv)
498	ENDIF
499	CALL xios_set_field_attr("other_field3D_resend_W",freq_op=freq_op_recv,freq_offset=freq_offset_recv)
500	ENDIF
501
502
503	CALL init_field2D(comm, other_params, lon, lat, domain_mask, field2D_init, &
504	x_lon, x_lat, x_mask, fieldX_init, y_lon, y_lat, y_mask, fieldY_init, fieldXY_init)
505
506	xy=size(domain_index)
507	x=size(X_index)
508	y=size(Y_index)
509	z=size(axis_index)
510	w=size(ensemble_value)
511
512	ALLOCATE(other_field2D(0:xy-1))
513	ALLOCATE(other_field3D(0:xy-1,0:z-1))
514	ALLOCATE(other_pressure(0:xy-1,0:z-1))
515	ALLOCATE(other_field3D_recv(0:xy-1,0:z-1))
516	ALLOCATE(other_field_Z(0:z-1))
517	ALLOCATE(other_field_X(0:x-1))
518	ALLOCATE(other_field_Y(0:y-1))
519	ALLOCATE(other_field_XY(0:x-1,0:y-1))
520	ALLOCATE(other_field_XYZ(0:x-1,0:y-1,0:z-1))
521	ALLOCATE(other_field_XZ(0:x-1,0:z-1))
522	ALLOCATE(other_field_YZ(0:y-1,0:z-1))
523
524	ALLOCATE(other_field2D_W(0:xy-1,0:w-1))
525	ALLOCATE(other_field3D_W(0:xy-1,0:z-1,0:w-1))
526	ALLOCATE(other_pressure_W(0:xy-1,0:z-1,0:w-1))
527	ALLOCATE(other_field3D_recv_W(0:xy-1,0:z-1,0:w-1))
528	ALLOCATE(other_field_ZW(0:z-1,0:w-1))
529	ALLOCATE(other_field_XW(0:x-1,0:w-1))
530	ALLOCATE(other_field_YW(0:y-1,0:w-1))
531	ALLOCATE(other_field_XYW(0:x-1,0:y-1,0:w-1))
532	ALLOCATE(other_field_XYZW(0:x-1,0:y-1,0:z-1,0:w-1))
533	ALLOCATE(other_field_XZW(0:x-1,0:z-1,0:w-1))
534	ALLOCATE(other_field_YZW(0:y-1,0:z-1,0:w-1))
535	ALLOCATE(other_field0D_W(0:w-1))
536
537
538
539	DO i=0,xy-1
540	IF (domain_index(i)/=-1) THEN
541	other_field2D(i)=field2D_init(domain_index(i))
542	ENDIF
543	ENDDO
544
545	DO i=0,x-1
546	IF (X_index(i)/=-1) THEN
547	other_field_X(i)=fieldX_init(X_index(i))
548	ENDIF
549	ENDDO
550
551	DO i=0,y-1
552	IF (Y_index(i)/=-1) THEN
553	other_field_Y(i)=fieldY_init(Y_index(i))
554	ENDIF
555	ENDDO
556
557	DO j=0,y-1
558	DO i=0,x-1
559	IF (Y_index(j)/=-1 .AND. X_index(i)/=-1) THEN
560	other_field_XY(i,j)=fieldXY_init(X_index(i),Y_index(j))
561	ENDIF
562	ENDDO
563	ENDDO
564
565	DO i=0,z-1
566	IF (axis_index(i)/=-1) THEN
567	k=axis_index(i)
568	IF (axis_mask(k)) THEN
569	other_field_Z(i)=axis_value(axis_index(k))
570	other_field3D(:,i)=other_field2D(:)*axis_value(axis_index(k))
571	other_field_XYZ(:,:,i)=other_field_XY(:,:)*axis_value(axis_index(k))
572	other_field_XZ(:,i)=other_field_X(:)*axis_value(axis_index(k))
573	other_field_YZ(:,i)=other_field_Y(:)*axis_value(axis_index(k))
574	ENDIF
575	ENDIF
576	ENDDO
577
578	pressure=1e20
579	DO j=0,z-1
580	pressure(:,j)=axis_value(j) ;
581	DO i=0,xy-1
582	IF (domain_index(i)/=-1) THEN
583	k=domain_index(i)
584	IF (domain_mask(k)) THEN
585	dist=sqrt((lat(k)/90.)2+(lon(k)/180.)2) ;
586	other_pressure(i,j)=other_pressure(i,j)(1+other_params%pressure_factorexp(-(4dist)*2))
587	ENDIF
588	ENDIF
589	ENDDO
590	ENDDO
591
592	other_field0D = 1
593
594
595	other_field2D_W(:,0) = other_field2D(:)(1+0.1ensemble_value(0))
596	other_field3D_W(:,:,0)= other_field3D(:,:)(1+0.1ensemble_value(0))
597	other_pressure_W(:,:,0) = other_pressure(:,:)(1+0.1ensemble_value(0))
598	other_field_XW(:,0) = other_field_X(:)(1+0.1ensemble_value(0))
599	other_field_YW(:,0) = other_field_Y(:)(1+0.1ensemble_value(0))
600	other_field_XYW(:,:,0) = other_field_XY(:,:)(1+0.1ensemble_value(0))
601	other_field_ZW(:,0) = other_field_Z(:)(1+0.1ensemble_value(0))
602	other_field_XYZW(:,:,:,0) = other_field_XYZ(:,:,:)(1+0.1ensemble_value(0))
603	other_field_XZW(:,:,0) = other_field_XZ(:,:)(1+0.1ensemble_value(0))
604	other_field_YZW(:,:,0) = other_field_YZ(:,:)(1+0.1ensemble_value(0))
605	other_field0D_W(0) = other_field0D(1+0.1ensemble_value(0))
606
607
608	ok_other_field2D=xios_is_valid_field("other_field2D") ;
609	ok_other_field3D=xios_is_valid_field("other_field3D") ;
610	ok_other_pressure=xios_is_valid_field("other_pressure") ;
611	ok_other_field2D_sub=xios_is_valid_field("other_field2D_sub") ;
612	ok_other_field3D_sub=xios_is_valid_field("other_field3D_sub") ;
613	ok_other_field_X=xios_is_valid_field("other_field_X") ;
614	ok_other_field_Y=xios_is_valid_field("other_field_Y") ;
615	ok_other_field_XY=xios_is_valid_field("other_field_XY") ;
616	ok_other_field_Z=xios_is_valid_field("other_field_Z") ;
617	ok_other_field_XYZ=xios_is_valid_field("other_field_XYZ") ;
618	ok_other_field_XZ=xios_is_valid_field("other_field_XZ") ;
619	ok_other_field_YZ=xios_is_valid_field("other_field_YZ") ;
620	ok_other_field0D=xios_is_valid_field("other_field0D") ;
621
622	ok_other_field2D_W=xios_is_valid_field("other_field2D_W") ;
623	ok_other_field3D_W=xios_is_valid_field("other_field3D_W") ;
624	ok_other_pressure_W=xios_is_valid_field("other_pressure_W") ;
625	ok_other_field2D_sub_W=xios_is_valid_field("other_field2D_sub_W") ;
626	ok_other_field3D_sub_W=xios_is_valid_field("other_field3D_sub_W") ;
627	ok_other_field_XW=xios_is_valid_field("other_field_XW") ;
628	ok_other_field_YW=xios_is_valid_field("other_field_YW") ;
629	ok_other_field_XYZW=xios_is_valid_field("other_field_XYW") ;
630	ok_other_field_ZW=xios_is_valid_field("other_field_ZW") ;
631	ok_other_field_XYZW=xios_is_valid_field("other_field_XYZW") ;
632	ok_other_field_XZW=xios_is_valid_field("other_field_XZW") ;
633	ok_other_field_YZW=xios_is_valid_field("other_field_YZW") ;
634	ok_other_field0D_W=xios_is_valid_field("other_field0D_W") ;
635
636
637	!!!!!!!!!!!!!!!!!!!!! end of duplicated section !!!!!!!!!!!!!!!!!!!!!!!!!!
638
639	ts=1
640	cdate=cdate+dtime
641	CALL xios_close_context_definition()
642	CALL xios_set_current_context(trim(model_id))
643
644
645	DO while ( cdate <= edate )
646
647	!!! Mise a jour du pas de temps
648	CALL xios_update_calendar(ts)
649
650	IF (ok_field2D) CALL xios_send_field("field2D",field2D)
651	IF (ok_field3D) CALL xios_send_field("field3D",field3D)
652	IF (ok_pressure) CALL xios_send_field("pressure",pressure)
653	IF (ok_field_X) CALL xios_send_field("field_X",field_X)
654	IF (ok_field_Y) CALL xios_send_field("field_Y",field_Y)
655	IF (ok_field_XY) CALL xios_send_field("field_XY",field_XY)
656	IF (ok_field_Z) CALL xios_send_field("field_Z",field_Z)
657	IF (ok_field_XYZ) CALL xios_send_field("field_XYZ",field_XYZ)
658	IF (ok_field_XZ) CALL xios_send_field("field_XZ",field_XZ)
659	IF (ok_field_YZ) CALL xios_send_field("field_YZ",field_YZ)
660	IF (ok_field0D) CALL xios_send_field("field0D",field0D)
661
662	IF ( MOD(params%field_sub_offset+ts-1,params%field_sub_freq)==0) THEN
663	IF (ok_field2D_sub) CALL xios_send_field("field2D_sub",field2D*10)
664	IF (ok_field3D_sub) CALL xios_send_field("field3D_sub",field3D*10)
665	ENDIF
666
667	IF (ok_field3D_recv) THEN
668	IF (xios_field_is_active("field3D_resend",.TRUE.)) THEN
669	CALL xios_recv_field("field3D_recv",field3D_recv)
670	CALL xios_send_field("field3D_resend",field3D_recv)
671	ENDIF
672	ENDIF
673
674	! ensemble
675	IF (ok_field2D_W) CALL xios_send_field("field2D_W",field2D_W)
676	IF (ok_field3D_W) CALL xios_send_field("field3D_W",field3D_W)
677	IF (ok_pressure_W) CALL xios_send_field("pressure_W",pressure_W)
678	IF (ok_field_XW) CALL xios_send_field("field_XW",field_XW)
679	IF (ok_field_YW) CALL xios_send_field("field_YW",field_YW)
680	IF (ok_field_XYW) CALL xios_send_field("field_XYW",field_XYW)
681	IF (ok_field_ZW) CALL xios_send_field("field_ZW",field_ZW)
682	IF (ok_field_XYZW) CALL xios_send_field("field_XYZW",field_XYZW)
683	IF (ok_field_XZW) CALL xios_send_field("field_XZW",field_XZW)
684	IF (ok_field_YZW) CALL xios_send_field("field_YZW",field_YZW)
685	IF (ok_field0D_W) CALL xios_send_field("field0D_W",field0D_W)
686
687	IF ( MOD(params%field_sub_offset+ts-1,params%field_sub_freq)==0) THEN
688	IF (ok_field2D_sub_W) CALL xios_send_field("field2D_sub_W",field2D_W*10)
689	IF (ok_field3D_sub_W) CALL xios_send_field("field3D_sub_W",field3D_W*10)
690	ENDIF
691
692	IF (ok_field3D_recv_W) THEN
693	IF (xios_field_is_active("field3D_resend_W",.TRUE.)) THEN
694	CALL xios_recv_field("field3D_recv_W",field3D_recv_W)
695	CALL xios_send_field("field3D_resend_W",field3D_recv_W)
696	ENDIF
697	ENDIF
698
699	field2D=field2D+1
700	field3D=field3D+1
701	field0D=field0D+1
702
703	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
704	! duplicated section for other_ !
705	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
706
707
708	IF (ok_other_field2D) CALL xios_send_field("other_field2D", other_field2D)
709	IF (ok_other_field3D) CALL xios_send_field("other_field3D", other_field3D)
710	IF (ok_other_pressure) CALL xios_send_field("other_pressure", other_pressure)
711	IF (ok_other_field_X) CALL xios_send_field("other_field_X", other_field_X)
712	IF (ok_other_field_Y) CALL xios_send_field("other_field_Y", other_field_Y)
713	IF (ok_other_field_XY) CALL xios_send_field("other_field_XY", other_field_XY)
714	IF (ok_other_field_Z) CALL xios_send_field("other_field_Z", other_field_Z)
715	IF (ok_other_field_XY) CALL xios_send_field("other_field_XYZ", other_field_XYZ)
716	IF (ok_other_field_XY) CALL xios_send_field("other_field_XZ", other_field_XZ)
717	IF (ok_other_field_XY) CALL xios_send_field("other_field_YZ", other_field_YZ)
718	IF (ok_other_field0D) CALL xios_send_field("other_field0D", other_field0D)
719
720	IF ( MOD(other_params%field_sub_offset+ts-1,other_params%field_sub_freq)==0) THEN
721	IF (ok_other_field2D_sub) CALL xios_send_field("other_field2D_sub",other_field2D*10)
722	IF (ok_other_field3D_sub) CALL xios_send_field("other_field3D_sub",other_field3D*10)
723	ENDIF
724
725	IF (ok_other_field3D_recv) THEN
726	IF (xios_field_is_active("other_field3D_resend",.TRUE.)) THEN
727	CALL xios_recv_field("other_field3D_recv",other_field3D_recv)
728	CALL xios_send_field("other_field3D_resend",other_field3D_recv)
729	ENDIF
730	ENDIF
731
732	! ensemble
733	IF (ok_other_field2D_W) CALL xios_send_field("other_field2D_W",other_field2D_W)
734	IF (ok_other_field3D_W) CALL xios_send_field("other_field3D_W",other_field3D_W)
735	IF (ok_other_pressure_W) CALL xios_send_field("other_pressure_W",other_pressure_W)
736	IF (ok_other_field_XW) CALL xios_send_field("other_field_XW",other_field_XW)
737	IF (ok_other_field_YW) CALL xios_send_field("other_field_YW",other_field_YW)
738	IF (ok_other_field_XYW) CALL xios_send_field("other_field_XYW",other_field_XYW)
739	IF (ok_other_field_ZW) CALL xios_send_field("other_field_ZW",other_field_ZW)
740	IF (ok_other_field_XYW) CALL xios_send_field("other_field_XYZW",other_field_XYZW)
741	IF (ok_other_field_XYW) CALL xios_send_field("other_field_XZW",other_field_XZW)
742	IF (ok_other_field_XYW) CALL xios_send_field("other_field_YZW",other_field_YZW)
743	IF (ok_other_field0D_W) CALL xios_send_field("other_field0D_W",other_field0D_W)
744
745	IF ( MOD(other_params%field_sub_offset+ts-1,other_params%field_sub_freq)==0) THEN
746	IF (ok_other_field2D_sub_W) CALL xios_send_field("other_field2D_sub_W",other_field2D_W*10)
747	IF (ok_other_field3D_sub_W) CALL xios_send_field("other_field3D_sub_W",other_field3D_W*10)
748	ENDIF
749
750	IF (ok_other_field3D_recv_W) THEN
751	IF (xios_field_is_active("other_field3D_resend_W",.TRUE.)) THEN
752	CALL xios_recv_field("other_field3D_recv_W",other_field3D_recv_W)
753	CALL xios_send_field("other_field3D_resend_W",other_field3D_recv_W)
754	ENDIF
755	ENDIF
756
757	other_field2D=other_field2D+1
758	other_field3D=other_field3D+1
759	other_field0D=other_field0D+1
760
761
762	!!!!!!!!!!!!!!!!!!!!!! end of duplicated section !!!!!!!!!!!!!!!!!
763
764
765
766	!! Boucle d'attente
767	CALL wait_us(int(elapsed_per_timestep*1.e6)) ! micro-secondes
768	cdate=cdate+dtime
769	ts=ts+1
770	ENDDO
771
772	CALL xios_context_finalize()
773	CALL MPI_COMM_FREE(comm, ierr)
774
775	END SUBROUTINE model
776
777	SUBROUTINE init_model_params(prefix,params)
778	IMPLICIT NONE
779	CHARACTER(LEN=*) :: prefix
780	TYPE(tmodel_params) :: params
781
782
783	IF (.NOT. xios_getvar(prefix//"timestep", params%timestep) ) params%timestep="1h"
784	IF (.NOT. xios_getvar(prefix//"domain", params%domain) ) params%domain="lmdz"
785	IF (.NOT. xios_getvar(prefix//"domain_mask", params%domain_mask) ) params%domain_mask=.FALSE.
786	IF (.NOT. xios_getvar(prefix//"domain_mask_type", params%domain_mask_type) ) params%domain_mask_type="cross"
787	IF (.NOT. xios_getvar(prefix//"scalar_mask", params%scalar_mask) ) params%scalar_mask=.FALSE.
788	IF (.NOT. xios_getvar(prefix//"scalar_mask_type", params%scalar_mask_type) ) params%scalar_mask_type="none"
789	IF (.NOT. xios_getvar(prefix//"axis", params%axis) ) params%axis="pressure"
790	IF (.NOT. xios_getvar(prefix//"axis_mask", params%axis_mask) ) params%axis_mask=.FALSE.
791	IF (.NOT. xios_getvar(prefix//"ni", params%ni) ) params%ni=36
792	IF (.NOT. xios_getvar(prefix//"nj", params%nj) ) params%nj=18
793	IF (.NOT. xios_getvar(prefix//"nlev", params%nlev) ) params%nlev=10
794	IF (.NOT. xios_getvar(prefix//"init_field2D", params%init_field2D) ) params%init_field2D="academic"
795	IF (.NOT. xios_getvar(prefix//"pressure_factor", params%pressure_factor) ) params%pressure_factor=0.
796	IF (.NOT. xios_getvar(prefix//"mask3d", params%mask3d) ) params%mask3d=.FALSE.
797	IF (.NOT. xios_getvar(prefix//"field_sub_freq", params%field_sub_freq) ) params%field_sub_freq=1
798	IF (.NOT. xios_getvar(prefix//"field_sub_offset", params%field_sub_offset) ) params%field_sub_offset=0
799
800	IF (.NOT. xios_getvar(prefix//"domain_proc_n", params%domain_proc_n)) params%domain_proc_n=0
801	IF (.NOT. xios_getvar(prefix//"axis_proc_n", params%axis_proc_n)) params%axis_proc_n=0
802	IF (.NOT. xios_getvar(prefix//"ensemble_proc_n", params%ensemble_proc_n)) params%ensemble_proc_n=1
803	IF ((.NOT. xios_getvar(prefix//"domain_proc_frac", params%domain_proc_frac)) .AND. &
804	(.NOT. xios_getvar("prefix//axis_proc_frac", params%axis_proc_frac))) THEN
805	params%domain_proc_frac=1.0
806	params%axis_proc_frac=0.0
807	ELSE IF (.NOT. xios_getvar(prefix//"domain_proc_frac", params%domain_proc_frac)) THEN
808	params%domain_proc_frac=0.0
809	ELSE IF (.NOT. xios_getvar(prefix//"axis_proc_frac", params%axis_proc_frac)) THEN
810	params%axis_proc_frac=0.0
811	ENDIF
812
813	END SUBROUTINE init_model_params
814
815
816	SUBROUTINE init_domain(domain_id, comm, params, return_ni, return_nj, &
817	return_lon, return_lat, return_mask, return_index, &
818	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
819	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
820	IMPLICIT NONE
821	CHARACTER(LEN=*),INTENT(IN) :: domain_id
822	TYPE(tmodel_params),INTENT(IN) :: params
823	INTEGER,INTENT(IN) :: comm
824	INTEGER,INTENT(OUT) :: return_ni
825	INTEGER,INTENT(OUT) :: return_nj
826	DOUBLE PRECISION, POINTER :: return_lon(:)
827	DOUBLE PRECISION, POINTER :: return_lat(:)
828	LOGICAL, POINTER :: return_mask(:)
829	INTEGER, POINTER :: return_index(:)
830	DOUBLE PRECISION, POINTER :: return_X_lon(:)
831	DOUBLE PRECISION, POINTER :: return_X_lat(:)
832	LOGICAL, POINTER :: return_X_mask(:)
833	INTEGER, POINTER :: return_X_index(:)
834	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
835	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
836	LOGICAL, POINTER :: return_Y_mask(:)
837	INTEGER, POINTER :: return_Y_index(:)
838
839	IF (params%domain=="lmdz") THEN
840	CALL init_domain_lmdz(domain_id,comm,params, return_ni, return_nj, &
841	return_lon, return_lat, return_mask, return_index, &
842	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
843	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
844
845	ELSE IF (params%domain=="orchidee") THEN
846	CALL init_domain_orchidee(domain_id,comm,params, return_ni, return_nj, &
847	return_lon, return_lat, return_mask, return_index, &
848	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
849	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
850
851	ELSE IF (params%domain=="nemo") THEN
852	CALL init_domain_nemo(domain_id,comm,params, return_ni, return_nj, &
853	return_lon, return_lat, return_mask, return_index, &
854	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
855	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
856
857	ELSE IF (params%domain=="dynamico") THEN
858	CALL init_domain_dynamico(domain_id,comm,params, return_ni, return_nj, &
859	return_lon, return_lat, return_mask, return_index, &
860	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
861	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
862
863	ELSE IF (params%domain=="gaussian") THEN
864	CALL init_domain_gaussian(domain_id,comm,params, return_ni, return_nj, &
865	return_lon, return_lat, return_mask, return_index, &
866	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
867	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
868	ENDIF
869
870	END SUBROUTINE init_domain
871
872	SUBROUTINE init_field2D(comm,params, lon, lat, mask, return_field, &
873	x_lon, x_lat, x_mask, return_fieldx, &
874	y_lon, y_lat, y_mask, return_fieldy, return_fieldXY)
875	IMPLICIT NONE
876	TYPE(tmodel_params) :: params
877	TYPE(xios_context) :: ctx_hdl
878	INTEGER :: comm
879	DOUBLE PRECISION, POINTER :: lon(:)
880	DOUBLE PRECISION, POINTER :: lat(:)
881	LOGICAL, POINTER :: mask(:)
882	DOUBLE PRECISION, POINTER :: return_field(:)
883
884	DOUBLE PRECISION, POINTER :: X_lon(:)
885	DOUBLE PRECISION, POINTER :: X_lat(:)
886	LOGICAL, POINTER :: X_mask(:)
887	DOUBLE PRECISION, POINTER :: return_fieldX(:)
888
889	DOUBLE PRECISION, POINTER :: Y_lon(:)
890	DOUBLE PRECISION, POINTER :: Y_lat(:)
891	LOGICAL, POINTER :: Y_mask(:)
892	DOUBLE PRECISION, POINTER :: return_fieldY(:)
893	DOUBLE PRECISION, POINTER :: return_fieldXY(:,:)
894
895	IF (params%init_field2D=="academic") THEN
896	CALL init_field2D_academic(comm,params, lon, lat, mask, return_field, X_lon, X_lat, X_mask, return_fieldX, &
897	Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
898	ELSE IF (params%init_field2D=="constant") THEN
899	CALL init_field2D_constant(comm,params, lon, lat, mask, return_field, X_lon, X_lat, X_mask, return_fieldX, &
900	Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
901	ELSE IF (params%init_field2D=="rank") THEN
902	CALL init_field2D_rank(comm,params, lon, lat, mask, return_field, X_lon, X_lat, X_mask, &
903	return_fieldX, Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
904	ENDIF
905
906	END SUBROUTINE init_field2D
907
908
909	SUBROUTINE init_axis(axis_id,comm,params, return_value, return_mask, return_index)
910	IMPLICIT NONE
911	CHARACTER(LEN=*) :: axis_id
912	TYPE(tmodel_params) :: params
913	TYPE(xios_context) :: ctx_hdl
914	INTEGER :: comm
915	DOUBLE PRECISION, POINTER :: return_value(:)
916	LOGICAL, POINTER :: return_mask(:)
917	INTEGER, POINTER :: return_index(:)
918
919	IF (params%axis=="pressure") THEN
920	CALL init_axis_pressure(axis_id,comm,params, return_value, return_mask, return_index)
921	ENDIF
922
923	END SUBROUTINE init_axis
924
925	SUBROUTINE init_ensemble(ensemble_id,comm,params, return_value)
926	IMPLICIT NONE
927	CHARACTER(LEN=*) :: ensemble_id
928	TYPE(tmodel_params) :: params
929	TYPE(xios_context) :: ctx_hdl
930	INTEGER :: comm
931	DOUBLE PRECISION, POINTER :: return_value(:)
932	INTEGER :: domain_proc_rank, domain_proc_size
933	INTEGER :: axis_proc_rank, axis_proc_size
934	INTEGER :: ensemble_proc_size, ensemble_proc_rank
935
936	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
937	ALLOCATE(return_value(0:0))
938	return_value(0)=ensemble_proc_rank
939
940	IF (xios_is_valid_axis(TRIM(ensemble_id))) THEN
941	CALL xios_set_axis_attr(ensemble_id, n_glo=ensemble_proc_size, begin=ensemble_proc_rank, n=1, value=return_value(:), unit='-')
942	ENDIF
943
944	END SUBROUTINE init_ensemble
945
946
947	SUBROUTINE init_domain_gaussian(domain_id, comm, params, return_ni, return_nj, &
948	return_lon,return_lat,return_mask, return_index, &
949	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
950	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
951	IMPLICIT NONE
952	CHARACTER(LEN=*) :: domain_id
953	TYPE(tmodel_params) :: params
954	TYPE(xios_context) :: ctx_hdl
955	INTEGER :: comm
956	INTEGER :: return_ni
957	INTEGER :: return_nj
958	DOUBLE PRECISION, POINTER :: return_lon(:)
959	DOUBLE PRECISION, POINTER :: return_lat(:)
960	LOGICAL, POINTER :: return_mask(:)
961	INTEGER, POINTER :: return_index(:)
962	DOUBLE PRECISION, POINTER :: return_X_lon(:)
963	DOUBLE PRECISION, POINTER :: return_X_lat(:)
964	LOGICAL, POINTER :: return_X_mask(:)
965	INTEGER, POINTER :: return_X_index(:)
966	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
967	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
968	LOGICAL, POINTER :: return_Y_mask(:)
969	INTEGER, POINTER :: return_Y_index(:)
970	INTEGER :: domain_proc_rank, domain_proc_size
971	INTEGER :: axis_proc_rank, axis_proc_size
972	INTEGER :: ensemble_proc_size, ensemble_proc_rank
973
974	INTEGER :: mpi_rank, mpi_size
975	INTEGER :: ierr
976
977	! INTEGER, PARAMETER :: nlon=60
978	! INTEGER, PARAMETER :: nlat=30
979	! INTEGER,PARAMETER :: ni_glo=100
980	! INTEGER,PARAMETER :: nj_glo=100
981	! INTEGER,PARAMETER :: llm=5
982	! DOUBLE PRECISION :: lval(llm)=1
983	! TYPE(xios_field) :: field_hdl
984	! TYPE(xios_fieldgroup) :: fieldgroup_hdl
985	! TYPE(xios_file) :: file_hdl
986	! LOGICAL :: ok
987
988	! DOUBLE PRECISION,ALLOCATABLE :: lon_glo(:),lat_glo(:)
989	! DOUBLE PRECISION,ALLOCATABLE :: bounds_lon_glo(:,:),bounds_lat_glo(:,:)
990	! DOUBLE PRECISION,ALLOCATABLE :: field_A_glo(:,:)
991	! INTEGER,ALLOCATABLE :: i_index_glo(:)
992	! INTEGER,ALLOCATABLE :: i_index(:)
993	! LOGICAL,ALLOCATABLE :: mask_glo(:),mask(:)
994	! INTEGER,ALLOCATABLE :: n_local(:),local_neighbor(:,:)
995	! DOUBLE PRECISION,ALLOCATABLE :: lon(:),lat(:),field_A_srf(:,:), lonvalue(:) ;
996	! DOUBLE PRECISION,ALLOCATABLE :: bounds_lon(:,:),bounds_lat(:,:) ;
997	! INTEGER :: ni,ibegin,iend,nj,jbegin,jend
998	! INTEGER :: i,j,l,ts,n, nbMax
999	! INTEGER :: ncell_glo,ncell,ind
1000	! REAL :: ilon,ilat
1001	! DOUBLE PRECISION, PARAMETER :: Pi=3.14159265359
1002	! INTEGER :: list_ind(nlon,nlat)
1003	! INTEGER :: rank,j1,j2,np,ncell_x
1004	! INTEGER :: data_n_index
1005	! INTEGER,ALLOCATABLE :: data_i_index(:)
1006	! DOUBLE PRECISION,ALLOCATABLE :: field_A_compressed(:,:)
1007
1008
1009	INTEGER :: nlon, nlat
1010	DOUBLE PRECISION,ALLOCATABLE :: lon_glo(:),lat_glo(:)
1011	DOUBLE PRECISION,ALLOCATABLE :: bounds_lon_glo(:,:),bounds_lat_glo(:,:)
1012	INTEGER,ALLOCATABLE :: i_index_glo(:)
1013	INTEGER,ALLOCATABLE :: i_index(:)
1014	LOGICAL,ALLOCATABLE :: mask_glo(:),mask(:)
1015	DOUBLE PRECISION,ALLOCATABLE :: lon(:),lat(:),field_A_srf(:,:), lonvalue(:) ;
1016	DOUBLE PRECISION,ALLOCATABLE :: bounds_lon(:,:),bounds_lat(:,:) ;
1017
1018	INTEGER :: ni,ibegin,iend,nj,jbegin,jend
1019	INTEGER :: i,j,l,ts,n, nbMax
1020	INTEGER :: ncell_glo,ncell,ind
1021	REAL :: ilon,ilat
1022	DOUBLE PRECISION, PARAMETER :: Pi=3.14159265359
1023	INTEGER,ALLOCATABLE :: list_ind(:,:)
1024	INTEGER :: rank,j1,j2,np,ncell_x
1025	INTEGER :: data_n_index
1026	INTEGER,ALLOCATABLE :: data_i_index(:)
1027
1028
1029
1030	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1031	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1032
1033	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
1034
1035	mpi_rank=domain_proc_rank
1036	mpi_size=domain_proc_size
1037	nlon=params%ni
1038	nlat=params%nj
1039
1040
1041	ncell_glo=0
1042	DO j=1,nlat
1043	n = NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1044	IF (n<8) n=8
1045	ncell_glo=ncell_glo+n
1046	ENDDO
1047
1048	ALLOCATE(lon_glo(ncell_glo))
1049	ALLOCATE(lat_glo(ncell_glo))
1050	ALLOCATE(bounds_lon_glo(4,ncell_glo))
1051	ALLOCATE(bounds_lat_glo(4,ncell_glo))
1052	ALLOCATE(i_index_glo(ncell_glo))
1053	ALLOCATE(mask_glo(ncell_glo))
1054	ALLOCATE(list_ind(nlon,nlat))
1055
1056	ind=0
1057	DO j=1,nlat
1058	n = NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1059	if (j==1) PRINT*,"--- ",n
1060	if (j==nlat) PRINT*,"--- ",n
1061	IF (n<8) n=8
1062
1063	DO i=1,n
1064	ind=ind+1
1065	list_ind(i,j)=ind
1066	ilon=i-0.5
1067	ilat=j-0.5
1068
1069	lat_glo(ind)= 90-(ilat*180./nlat)
1070	lon_glo(ind)= (ilon*360./n)-180.
1071
1072
1073	bounds_lat_glo(1,ind)= 90-((ilat-0.5)*180./nlat)
1074	bounds_lon_glo(1,ind)=((ilon-0.5)*360./n)-180.
1075
1076	bounds_lat_glo(2,ind)= 90-((ilat-0.5)*180./nlat)
1077	bounds_lon_glo(2,ind)=((ilon+0.5)*360./n)-180.
1078
1079	bounds_lat_glo(3,ind)= 90-((ilat+0.5)*180./nlat)
1080	bounds_lon_glo(3,ind)=((ilon+0.5)*360./n)-180.
1081
1082	bounds_lat_glo(4,ind)= 90-((ilat+0.5)*180./nlat)
1083	bounds_lon_glo(4,ind)=((ilon-0.5)*360./n)-180.
1084
1085	ENDDO
1086	ENDDO
1087
1088	! mpi_size=32
1089	rank=(mpi_size-1)/2
1090	ncell_x=sqrt(ncell_glo*1./mpi_size)
1091
1092	j1=nlat/2
1093	DO WHILE(rank>=0)
1094	j2=MAX(j1-ncell_x+1,1)
1095	j=(j1+j2)/2
1096	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1097	np = MIN(n/ncell_x,rank+1) ;
1098	if (j2==1) np=rank+1
1099
1100	PRINT *,"domain ",j2,j1,rank,np ;
1101	DO j=j2,j1
1102	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1103	IF (n<8) n=8
1104	DO i=1,n
1105	ind=list_ind(i,j)
1106	IF ( (i-1) < MOD(n,np)*(n/np+1)) THEN
1107	i_index_glo(ind) = rank - (i-1)/(n/np+1)
1108	ELSE
1109	i_index_glo(ind) = rank-(MOD(n,np)+ (i-1-MOD(n,np)*(n/np+1))/(n/np))
1110	ENDIF
1111	ENDDO
1112	ENDDO
1113	rank=rank-np
1114	j1=j2-1
1115	ENDDO
1116
1117	rank=(mpi_size-1)/2+1
1118	ncell_x=sqrt(ncell_glo*1./mpi_size)
1119
1120	j1=nlat/2+1
1121	DO WHILE(rank<=mpi_size-1)
1122	j2=MIN(j1+ncell_x-1,nlat)
1123	j=(j1+j2)/2
1124	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1125	np = MIN(n/ncell_x,mpi_size-rank) ;
1126	if (j2==nlat) np=mpi_size-rank
1127
1128	PRINT *,"domain ",j2,j1,rank,np ;
1129	DO j=j1,j2
1130	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
1131	IF (n<8) n=8
1132	DO i=1,n
1133	ind=list_ind(i,j)
1134	IF ( (i-1) < MOD(n,np)*(n/np+1)) THEN
1135	i_index_glo(ind) = rank + (i-1)/(n/np+1)
1136	ELSE
1137	i_index_glo(ind) = rank+(MOD(n,np)+ (i-1-MOD(n,np)*(n/np+1))/(n/np))
1138	ENDIF
1139	ENDDO
1140	ENDDO
1141	rank=rank+np
1142	j1=j2+1
1143	ENDDO
1144
1145	ncell=0
1146	DO ind=1,ncell_glo
1147	IF (i_index_glo(ind)==mpi_rank) ncell=ncell+1
1148	ENDDO
1149	ALLOCATE(i_index(ncell))
1150	ALLOCATE(lon(ncell))
1151	ALLOCATE(lat(ncell))
1152	ALLOCATE(bounds_lon(4,ncell))
1153	ALLOCATE(bounds_lat(4,ncell))
1154	! ALLOCATE(field_A_srf(ncell,llm))
1155	ALLOCATE(mask(ncell))
1156	ncell=0
1157	data_n_index=0
1158	DO ind=1,ncell_glo
1159	IF (i_index_glo(ind)==mpi_rank) THEN
1160	ncell=ncell+1
1161	i_index(ncell)=ind-1
1162	lon(ncell)=lon_glo(ind)
1163	lat(ncell)=lat_glo(ind)
1164	bounds_lon(:,ncell)=bounds_lon_glo(:,ind)
1165	bounds_lat(:,ncell)=bounds_lat_glo(:,ind)
1166	field_A_srf(ncell,:)=i_index_glo(ind)
1167	IF (MOD(ind,8)>=0 .AND. MOD(ind,8)<2) THEN
1168	! mask(ncell)=.FALSE.
1169
1170	mask(ncell)=.TRUE.
1171	data_n_index=data_n_index+1
1172	ELSE
1173	mask(ncell)=.TRUE.
1174	data_n_index=data_n_index+1
1175	ENDIF
1176	ENDIF
1177	ENDDO
1178
1179	! ALLOCATE(field_A_compressed(data_n_index,llm))
1180	ALLOCATE(data_i_index(data_n_index))
1181	data_n_index=0
1182	DO ind=1,ncell
1183	IF (mask(ind)) THEN
1184	data_n_index=data_n_index+1
1185	data_i_index(data_n_index)=ind-1
1186	! field_A_compressed(data_n_index,:)=field_A_srf(ind,:)
1187	ENDIF
1188	ENDDO
1189
1190	ALLOCATE(return_lon(0:ncell-1))
1191	ALLOCATE(return_lat(0:ncell-1))
1192	ALLOCATE(return_mask(0:ncell-1))
1193	ALLOCATE(return_index(0:data_n_index-1))
1194	return_lon(0:ncell-1)=lon(1:ncell)
1195	return_lat(0:ncell-1)=lat(1:ncell)
1196	return_index(0:data_n_index-1)= data_i_index(1:data_n_index)
1197	CALL set_domain_mask(params, lon, lat, return_mask)
1198
1199
1200	ALLOCATE(return_X_lon(0:ncell-1))
1201	ALLOCATE(return_X_lat(0:ncell-1))
1202	ALLOCATE(return_X_mask(0:ncell-1))
1203	ALLOCATE(return_X_index(0:ncell-1))
1204
1205	return_X_lon=return_lon
1206	return_X_lat=return_lat
1207	return_X_index=return_index
1208	CALL set_domain_mask(params, return_X_lon,return_X_lat, return_X_mask)
1209
1210
1211	ALLOCATE(return_Y_lon(0:0))
1212	ALLOCATE(return_Y_lat(0:0))
1213	ALLOCATE(return_Y_mask(0:0))
1214	ALLOCATE(return_Y_index(0:0))
1215
1216	return_Y_lon(0)=lon_glo(ncell_glo/2)
1217	return_Y_lat(0)=lat_glo(ncell_glo/2)
1218	CALL set_domain_mask(params, return_Y_lon,return_Y_lat, return_Y_mask)
1219	return_Y_index(0)=0
1220
1221	return_ni=ncell
1222	return_nj=1
1223
1224
1225	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1226	CALL xios_set_domain_attr(TRIM(domain_id), type="unstructured", ni_glo=ncell_glo, ni=ncell, ibegin=0, i_index=i_index)
1227	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=1, data_ni=data_n_index, data_i_index=data_i_index, mask_1d=return_mask)
1228	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_1D=lon, latvalue_1D=lat, nvertex=4, bounds_lon_1D=bounds_lon, bounds_lat_1D=bounds_lat)
1229	ENDIF
1230
1231	IF (xios_is_valid_axis(TRIM(domain_id)//"_X")) THEN
1232	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ncell_glo, n=ncell, begin=0, index=i_index,value=return_X_lon)
1233	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", data_n=data_n_index, data_index=data_i_index, mask=return_X_mask)
1234	ENDIF
1235
1236	IF (xios_is_valid_axis(TRIM(domain_id)//"_Y")) THEN
1237	CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=1, begin=0, n=1, value=return_Y_lat, mask=return_Y_mask)
1238	ENDIF
1239
1240
1241
1242	END SUBROUTINE init_domain_gaussian
1243
1244
1245
1246	SUBROUTINE init_domain_dynamico(domain_id, comm, params, return_ni, return_nj, &
1247	return_lon,return_lat,return_mask, return_index, &
1248	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
1249	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
1250	IMPLICIT NONE
1251	CHARACTER(LEN=*) :: domain_id
1252	TYPE(tmodel_params) :: params
1253	TYPE(xios_context) :: ctx_hdl
1254	INTEGER :: comm
1255	INTEGER :: return_ni
1256	INTEGER :: return_nj
1257	DOUBLE PRECISION, POINTER :: return_lon(:)
1258	DOUBLE PRECISION, POINTER :: return_lat(:)
1259	LOGICAL, POINTER :: return_mask(:)
1260	INTEGER, POINTER :: return_index(:)
1261	DOUBLE PRECISION, POINTER :: return_X_lon(:)
1262	DOUBLE PRECISION, POINTER :: return_X_lat(:)
1263	LOGICAL, POINTER :: return_X_mask(:)
1264	INTEGER, POINTER :: return_X_index(:)
1265	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
1266	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
1267	LOGICAL, POINTER :: return_Y_mask(:)
1268	INTEGER, POINTER :: return_Y_index(:)
1269	INTEGER :: domain_proc_rank, domain_proc_size
1270	INTEGER :: axis_proc_rank, axis_proc_size
1271	INTEGER :: ensemble_proc_size, ensemble_proc_rank
1272
1273	INTEGER :: mpi_rank, mpi_size
1274	INTEGER :: ierr
1275	INTEGER :: ni,nj,ni_glo,nj_glo,nvertex
1276	INTEGER :: ibegin,jbegin
1277	INTEGER :: nbp,nbp_glo, offset
1278	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), lon(:), lat(:)
1279	DOUBLE PRECISION, ALLOCATABLE :: bounds_lon_glo(:,:), bounds_lat_glo(:,:), bounds_lon(:,:), bounds_lat(:,:)
1280	LOGICAL,ALLOCATABLE :: mask(:)
1281	LOGICAL,ALLOCATABLE :: dom_mask(:)
1282	INTEGER :: i,j
1283	INTEGER,ALLOCATABLE :: ibegin_all(:), ni_all(:)
1284
1285	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1286	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1287
1288	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
1289
1290	CALL xios_get_current_context(ctx_hdl)
1291
1292	CALL xios_context_initialize("grid_dynamico",comm)
1293	CALL xios_close_context_definition()
1294	CALL xios_get_domain_attr("domain_dynamico",ni_glo=ni_glo,nj_glo=nj_glo,ni=ni,nj=nj, ibegin=ibegin, jbegin=jbegin ,nvertex=nvertex)
1295	ALLOCATE(lon(ni),lat(ni),bounds_lon(nvertex,ni),bounds_lat(nvertex,ni))
1296	CALL xios_get_domain_attr("domain_dynamico", lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
1297	CALL xios_context_finalize
1298
1299	CALL xios_set_current_context(ctx_hdl)
1300
1301
1302	ALLOCATE(ibegin_all(mpi_size))
1303	ALLOCATE(ni_all(mpi_size))
1304	ALLOCATE(lon_glo(0:ni_glo-1))
1305	ALLOCATE(lat_glo(0:ni_glo-1))
1306	ALLOCATE(bounds_lon_glo(nvertex,0:ni_glo-1))
1307	ALLOCATE(bounds_lat_glo(nvertex,0:ni_glo-1))
1308
1309	CALL MPI_Allgather(ibegin, 1, MPI_INTEGER, ibegin_all, 1, MPI_INTEGER, comm, ierr)
1310	CALL MPI_Allgather(ni, 1, MPI_INTEGER, ni_all, 1, MPI_INTEGER, comm, ierr)
1311
1312	CALL MPI_AllgatherV(lon, ni, MPI_REAL8, lon_glo, ni_all, ibegin_all, MPI_REAL8, comm, ierr)
1313	CALL MPI_AllgatherV(lat, ni, MPI_REAL8, lat_glo, ni_all, ibegin_all, MPI_REAL8, comm, ierr)
1314	CALL MPI_AllgatherV(bounds_lon, ninvertex, MPI_REAL8, bounds_lon_glo, ni_allnvertex, ibegin_all*nvertex, MPI_REAL8, comm, ierr)
1315	CALL MPI_AllgatherV(bounds_lat, ninvertex, MPI_REAL8, bounds_lat_glo, ni_allnvertex, ibegin_all*nvertex, MPI_REAL8, comm, ierr)
1316
1317
1318	nbp_glo=ni_glo
1319	nbp=nbp_glo/domain_proc_size
1320	IF (domain_proc_rank < MOD(nbp_glo,domain_proc_size) ) THEN
1321	nbp=nbp+1
1322	offset=nbp*domain_proc_rank
1323	ELSE
1324	offset=nbp*domain_proc_rank + MOD(nbp_glo, domain_proc_size)
1325	ENDIF
1326
1327	DEALLOCATE(lat,lon,bounds_lon,bounds_lat)
1328	ALLOCATE(lat(0:nbp-1))
1329	ALLOCATE(lon(0:nbp-1))
1330	ALLOCATE(bounds_lon(nvertex,0:nbp-1))
1331	ALLOCATE(bounds_lat(nvertex,0:nbp-1))
1332	ALLOCATE(return_lon(0:nbp-1))
1333	ALLOCATE(return_lat(0:nbp-1))
1334	ALLOCATE(return_index(0:nbp-1))
1335
1336	DO i=0,nbp-1
1337	lat(i)=lat_glo(i+offset)
1338	lon(i)=lon_glo(i+offset)
1339	bounds_lon(:,i) = bounds_lon_glo(:,i+offset)
1340	bounds_lat(:,i) = bounds_lat_glo(:,i+offset)
1341	return_index(i)=i
1342	ENDDO
1343	return_lon=lon
1344	return_lat=lat
1345
1346	ALLOCATE(return_mask(0:nbp-1))
1347	CALL set_domain_mask(params, lon, lat, return_mask)
1348
1349
1350	ALLOCATE(return_X_lon(0:nbp-1))
1351	ALLOCATE(return_X_lat(0:nbp-1))
1352	ALLOCATE(return_X_mask(0:nbp-1))
1353	ALLOCATE(return_X_index(0:nbp-1))
1354
1355	DO i=0,nbp-1
1356	return_X_lon(i)=lon_glo(i+offset)
1357	return_X_lat(i)=lat_glo(i+offset)
1358	ENDDO
1359	CALL set_domain_mask(params, return_X_lon,return_X_lat, return_X_mask)
1360	DO i=0,nbp-1
1361	return_X_index(i)=i
1362	ENDDO
1363
1364	ALLOCATE(return_Y_lon(0:0))
1365	ALLOCATE(return_Y_lat(0:0))
1366	ALLOCATE(return_Y_mask(0:0))
1367	ALLOCATE(return_Y_index(0:0))
1368
1369	return_Y_lon(0)=lon_glo(nbp_glo/2)
1370	return_Y_lat(0)=lat_glo(nbp_glo/2)
1371	CALL set_domain_mask(params, return_Y_lon,return_Y_lat, return_Y_mask)
1372	return_Y_index(0)=0
1373
1374
1375
1376	return_ni=nbp
1377	return_nj=1
1378
1379	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1380	CALL xios_set_domain_attr(TRIM(domain_id), type="unstructured", ni_glo=ni_glo, ibegin=offset, ni=nbp, nvertex=nvertex)
1381	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=1, lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat, mask_1d=return_mask)
1382	ENDIF
1383
1384	IF (xios_is_valid_axis(TRIM(domain_id)//"_X")) THEN
1385	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ni_glo, begin=offset, n=nbp, value=return_X_lon)
1386	ENDIF
1387
1388	IF (xios_is_valid_axis(TRIM(domain_id)//"_Y")) THEN
1389	CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=1, begin=0, n=1, value=return_Y_lat, mask=return_Y_mask)
1390	ENDIF
1391
1392	END SUBROUTINE init_domain_dynamico
1393
1394
1395
1396
1397	SUBROUTINE init_domain_lmdz(domain_id, comm, params, return_ni, return_nj, &
1398	return_lon,return_lat,return_mask, return_index, &
1399	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
1400	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
1401	IMPLICIT NONE
1402	CHARACTER(LEN=*) :: domain_id
1403	TYPE(tmodel_params) :: params
1404	TYPE(xios_context) :: ctx_hdl
1405	INTEGER :: comm
1406	INTEGER :: return_ni
1407	INTEGER :: return_nj
1408	DOUBLE PRECISION, POINTER :: return_lon(:)
1409	DOUBLE PRECISION, POINTER :: return_lat(:)
1410	LOGICAL, POINTER :: return_mask(:)
1411	INTEGER, POINTER :: return_index(:)
1412	DOUBLE PRECISION, POINTER :: return_X_lon(:)
1413	DOUBLE PRECISION, POINTER :: return_X_lat(:)
1414	LOGICAL, POINTER :: return_X_mask(:)
1415	INTEGER, POINTER :: return_X_index(:)
1416	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
1417	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
1418	LOGICAL, POINTER :: return_Y_mask(:)
1419	INTEGER, POINTER :: return_Y_index(:)
1420	INTEGER :: domain_proc_rank, domain_proc_size
1421	INTEGER :: axis_proc_rank, axis_proc_size
1422	INTEGER :: ensemble_proc_size, ensemble_proc_rank
1423
1424	INTEGER :: mpi_rank, mpi_size
1425	INTEGER :: ierr
1426	INTEGER :: ni,nj,ni_glo,nj_glo
1427	INTEGER :: ibegin,jbegin
1428	INTEGER :: nbp,nbp_glo, offset
1429	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), lon(:), lat(:)
1430	LOGICAL,ALLOCATABLE :: mask(:)
1431	LOGICAL,ALLOCATABLE :: dom_mask(:)
1432	INTEGER :: i,j
1433
1434	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1435	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1436
1437	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
1438	ni_glo=params%ni
1439	nj_glo=params%nj
1440	nbp_glo=ni_glo*nj_glo
1441	nbp=nbp_glo/domain_proc_size
1442	IF (domain_proc_rank < MOD(nbp_glo,domain_proc_size) ) THEN
1443	nbp=nbp+1
1444	offset=nbp*domain_proc_rank
1445	ELSE
1446	offset=nbp*domain_proc_rank + MOD(nbp_glo, domain_proc_size)
1447	ENDIF
1448
1449
1450	ibegin=0 ; ni=ni_glo
1451	jbegin=offset / ni_glo
1452
1453	nj = (offset+nbp-1) / ni_glo - jbegin + 1
1454
1455	offset=offset-jbegin*ni
1456
1457	ALLOCATE(lon(0:ni-1), lat(0:nj-1), mask(0:ninj-1), dom_mask(0:ninj-1))
1458	mask(:)=.FALSE.
1459	mask(offset:offset+nbp-1)=.TRUE.
1460
1461	ALLOCATE(lon_glo(0:ni_glo-1), lat_glo(0:nj_glo-1))
1462
1463	DO i=0,ni_glo-1
1464	lon_glo(i)=-180+(i+0.5)*(360./ni_glo)
1465	ENDDO
1466
1467	DO j=0,nj_glo-1
1468	lat_glo(j)=-90+(j+0.5)*(180./nj_glo)
1469	ENDDO
1470
1471	lon(:)=lon_glo(:)
1472	lat(:)=lat_glo(jbegin:jbegin+nj-1)
1473
1474	ALLOCATE(return_lon(0:ni*nj-1))
1475	ALLOCATE(return_lat(0:ni*nj-1))
1476	ALLOCATE(return_mask(0:ni*nj-1))
1477	ALLOCATE(return_index(0:ni*nj-1))
1478
1479	ALLOCATE(return_X_lon(0:ni-1))
1480	ALLOCATE(return_X_lat(0:ni-1))
1481	ALLOCATE(return_X_mask(0:ni-1))
1482	ALLOCATE(return_X_index(0:ni-1))
1483
1484	ALLOCATE(return_Y_lon(0:nj-1))
1485	ALLOCATE(return_Y_lat(0:nj-1))
1486	ALLOCATE(return_Y_mask(0:nj-1))
1487	ALLOCATE(return_Y_index(0:nj-1))
1488
1489	DO j=0,nj-1
1490	DO i=0,ni-1
1491	return_lon(i+j*ni)=lon(i)
1492	return_lat(i+j*ni)=lat(j)
1493	return_index(i+jni)=i+jni
1494	ENDDO
1495	ENDDO
1496
1497	CALL set_domain_mask(params, return_lon,return_lat, dom_mask)
1498
1499	return_mask = mask .AND. dom_mask
1500
1501	return_X_lat(:)=lat_glo(nj_glo/2)
1502	return_X_lon(:)=lon_glo(:)
1503	CALL set_domain_mask(params, return_X_lon,return_X_lat, return_X_mask)
1504	DO i=0,ni-1
1505	return_X_index(i)=i
1506	ENDDO
1507
1508	return_Y_lon(:)=lon_glo(ni_glo/2)
1509	return_Y_lat(:)=lat_glo(jbegin:jbegin+nj-1)
1510	CALL set_domain_mask(params, return_Y_lon,return_Y_lat, return_Y_mask)
1511
1512	DO j=0,nj-1
1513	return_Y_index(j)=j
1514	ENDDO
1515
1516	return_ni=ni
1517	return_nj=nj
1518
1519	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1520	CALL xios_set_domain_attr(TRIM(domain_id), type="rectilinear", ni_glo=ni_glo, ibegin=ibegin, ni=ni, nj_glo=nj_glo, jbegin=jbegin, nj=nj)
1521	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=2, lonvalue_1d=lon, latvalue_1d=lat, mask_1d=return_mask)
1522	ENDIF
1523
1524	IF (xios_is_valid_axis(TRIM(domain_id)//"_X")) THEN
1525	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ni_glo, begin=ibegin, n=ni, value=return_X_lon)
1526	ENDIF
1527
1528	IF (xios_is_valid_axis(TRIM(domain_id)//"_Y")) THEN
1529	CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=nj_glo, begin=jbegin, n=nj, value=return_Y_lat)
1530	ENDIF
1531
1532	END SUBROUTINE init_domain_lmdz
1533
1534
1535	SUBROUTINE init_domain_orchidee(domain_id, comm, params, return_ni, return_nj, &
1536	return_lon,return_lat,return_mask, return_index, &
1537	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
1538	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
1539	IMPLICIT NONE
1540	CHARACTER(LEN=*) :: domain_id
1541	TYPE(tmodel_params) :: params
1542	TYPE(xios_context) :: ctx_hdl
1543	INTEGER :: comm
1544	INTEGER :: return_ni
1545	INTEGER :: return_nj
1546	DOUBLE PRECISION, POINTER :: return_lon(:)
1547	DOUBLE PRECISION, POINTER :: return_lat(:)
1548	LOGICAL, POINTER :: return_mask(:)
1549	INTEGER, POINTER :: return_index(:)
1550	DOUBLE PRECISION, POINTER :: return_X_lon(:)
1551	DOUBLE PRECISION, POINTER :: return_X_lat(:)
1552	LOGICAL, POINTER :: return_X_mask(:)
1553	INTEGER, POINTER :: return_X_index(:)
1554	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
1555	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
1556	LOGICAL, POINTER :: return_Y_mask(:)
1557	INTEGER, POINTER :: return_Y_index(:)
1558	INTEGER :: domain_proc_rank, domain_proc_size
1559	INTEGER :: axis_proc_rank, axis_proc_size
1560	INTEGER :: ensemble_proc_size, ensemble_proc_rank
1561
1562	INTEGER :: mpi_rank, mpi_size
1563	INTEGER :: ierr
1564	INTEGER :: ni,nj,ni_glo,nj_glo
1565	INTEGER :: ibegin,jbegin
1566	INTEGER :: nbp,nbp_glo, offset
1567	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), lon(:), lat(:)
1568	LOGICAL,ALLOCATABLE :: mask(:)
1569	INTEGER :: i,j, ij, pos, i_glo, j_glo, ij_begin, ij_end
1570
1571	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1572	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1573
1574	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
1575	ni_glo=params%ni
1576	nj_glo=params%nj
1577	nbp_glo=ni_glo*nj_glo
1578	nbp=nbp_glo/domain_proc_size
1579	IF (domain_proc_rank < MOD(nbp_glo,domain_proc_size) ) THEN
1580	nbp=nbp+1
1581	offset=nbp*domain_proc_rank
1582	ELSE
1583	offset=nbp*domain_proc_rank + MOD(nbp_glo, domain_proc_size)
1584	ENDIF
1585
1586
1587	ibegin=0 ; ni=ni_glo
1588	jbegin=offset / ni_glo
1589
1590	nj = (offset+nbp-1) / ni_glo - jbegin + 1
1591
1592	offset=offset-jbegin*ni
1593
1594	ij_begin=offset+jbegin*ni
1595	ij_end=ij_begin+nbp-1
1596
1597	ALLOCATE(lon(0:ni-1), lat(0:nj-1), mask(0:ni*nj-1))
1598	mask(:)=.FALSE.
1599	mask(offset:offset+nbp-1)=.TRUE.
1600
1601	ALLOCATE(lon_glo(0:ni_glo-1), lat_glo(0:nj_glo-1))
1602
1603	DO i=0,ni_glo-1
1604	lon_glo(i)=-180+(i+0.5)*(360./ni_glo)
1605	ENDDO
1606
1607	DO j=0,nj_glo-1
1608	lat_glo(j)=-90+(j+0.5)*(180./nj_glo)
1609	ENDDO
1610
1611	lon(:)=lon_glo(:)
1612	lat(:)=lat_glo(jbegin:jbegin+nj-1)
1613
1614	ALLOCATE(return_lon(0:ni*nj-1))
1615	ALLOCATE(return_lat(0:ni*nj-1))
1616	ALLOCATE(return_mask(0:ni*nj-1))
1617
1618	ALLOCATE(return_X_lon(0:ni-1))
1619	ALLOCATE(return_X_lat(0:ni-1))
1620	ALLOCATE(return_X_mask(0:ni-1))
1621
1622	ALLOCATE(return_Y_lon(0:nj-1))
1623	ALLOCATE(return_Y_lat(0:nj-1))
1624	ALLOCATE(return_Y_mask(0:nj-1))
1625
1626
1627	DO j=0,nj-1
1628	DO i=0,ni-1
1629	return_lon(i+j*ni)=lon(i)
1630	return_lat(i+j*ni)=lat(j)
1631	ENDDO
1632	ENDDO
1633
1634	pos=0
1635	DO j=0,nj-1
1636	DO i=0,ni-1
1637	i_glo=i
1638	j_glo=jbegin+j
1639	ij=j_glo*ni_glo+i_glo
1640	IF ( ij>=ij_begin .AND. ij<=ij_end) THEN
1641	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1642	pos=pos+1
1643	ENDIF
1644	ENDDO
1645	ENDDO
1646
1647	ALLOCATE(return_index(0:pos-1))
1648	return_index(:)=-1
1649	pos=0
1650	DO j=0,nj-1
1651	DO i=0,ni-1
1652	i_glo=i
1653	j_glo=jbegin+j
1654	ij=j_glo*ni_glo+i_glo
1655	IF ( ij>=ij_begin .AND. ij<=ij_end) THEN
1656	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1657	return_index(pos)=i+j*ni
1658	pos=pos+1
1659	ENDIF
1660	ENDDO
1661	ENDDO
1662
1663	CALL set_domain_mask(params, return_lon,return_lat, mask)
1664
1665	return_mask=mask
1666
1667	return_X_lat(:)=lat_glo(nj_glo/2)
1668	return_X_lon(:)=lon_glo(:)
1669	CALL set_domain_mask(params, return_X_lon,return_X_lat, return_X_mask)
1670
1671	pos=0
1672	DO i=0,ni-1
1673	i_glo=i
1674	j_glo=nj_glo/2
1675	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1676	pos=pos+1
1677	ENDDO
1678
1679	ALLOCATE(return_X_index(0:pos-1))
1680	return_X_index(:)=-1
1681	pos=0
1682	DO i=0,ni-1
1683	i_glo=i
1684	j_glo=nj_glo/2
1685	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1686	return_X_index(pos)=i
1687	pos=pos+1
1688	ENDDO
1689
1690
1691
1692	return_Y_lon(:)=lon_glo(ni_glo/2)
1693	return_Y_lat(:)=lat_glo(jbegin:jbegin+nj-1)
1694	CALL set_domain_mask(params, return_Y_lon,return_Y_lat, return_Y_mask)
1695
1696	pos=0
1697	DO j=0,nj-1
1698	i_glo=ni_glo/2
1699	j_glo=j+jbegin
1700	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1701	pos=pos+1
1702	ENDDO
1703
1704	ALLOCATE(return_Y_index(0:pos-1))
1705	return_Y_index=-1
1706	pos=0
1707	DO j=0,nj-1
1708	i_glo=ni_glo/2
1709	j_glo=j+jbegin
1710	IF ((MOD(j_glo,3)==1 .OR. MOD(j_glo,3)==2) .AND. (MOD(i_glo,5)==3 .OR. MOD(i_glo,5)==4)) CYCLE
1711	return_Y_index(pos)=j
1712	pos=pos+1
1713	ENDDO
1714
1715	return_ni=ni
1716	return_nj=nj
1717
1718	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1719	CALL xios_set_domain_attr(TRIM(domain_id), type="rectilinear", ni_glo=ni_glo, ibegin=ibegin, ni=ni, nj_glo=nj_glo, jbegin=jbegin, nj=nj)
1720	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=1, data_ni=size(return_index), data_i_index=return_index)
1721	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_1d=lon, latvalue_1d=lat, mask_1d=mask)
1722	ENDIF
1723
1724	IF (xios_is_valid_axis(TRIM(domain_id)//"_X")) THEN
1725	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ni_glo, begin=ibegin, n=ni, value=return_X_lon, data_n=size(return_X_index), data_index=return_X_index)
1726	ENDIF
1727
1728	IF (xios_is_valid_axis(TRIM(domain_id)//"_Y")) THEN
1729	CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=nj_glo, begin=jbegin, n=nj, value=return_Y_lat, data_n=size(return_Y_index), data_index=return_Y_index)
1730	ENDIF
1731
1732	END SUBROUTINE init_domain_orchidee
1733
1734
1735
1736	SUBROUTINE init_domain_nemo(domain_id, comm, params, return_ni, return_nj, &
1737	return_lon,return_lat,return_mask, return_index, &
1738	return_X_lon,return_X_lat, return_X_mask, return_X_index, &
1739	return_Y_lon,return_Y_lat, return_Y_mask, return_Y_index)
1740	IMPLICIT NONE
1741	CHARACTER(LEN=*) :: domain_id
1742	TYPE(tmodel_params) :: params
1743	TYPE(xios_context) :: ctx_hdl
1744	INTEGER :: comm
1745	INTEGER :: return_ni
1746	INTEGER :: return_nj
1747	DOUBLE PRECISION, POINTER :: return_lon(:)
1748	DOUBLE PRECISION, POINTER :: return_lat(:)
1749	LOGICAL, POINTER :: return_mask(:)
1750	INTEGER, POINTER :: return_index(:)
1751	DOUBLE PRECISION, POINTER :: return_X_lon(:)
1752	DOUBLE PRECISION, POINTER :: return_X_lat(:)
1753	LOGICAL, POINTER :: return_X_mask(:)
1754	INTEGER, POINTER :: return_X_index(:)
1755	DOUBLE PRECISION, POINTER :: return_Y_lon(:)
1756	DOUBLE PRECISION, POINTER :: return_Y_lat(:)
1757	LOGICAL, POINTER :: return_Y_mask(:)
1758	INTEGER, POINTER :: return_Y_index(:)
1759	INTEGER :: domain_proc_rank, domain_proc_size
1760	INTEGER :: axis_proc_rank, axis_proc_size
1761	INTEGER :: ensemble_proc_size, ensemble_proc_rank
1762
1763	INTEGER :: mpi_rank, mpi_size
1764	INTEGER :: ierr
1765	INTEGER :: ni,nj,ni_glo,nj_glo
1766	INTEGER :: ibegin,jbegin
1767	INTEGER :: offset_i, offset_j
1768	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), bounds_lon_glo(:,:), bounds_lat_glo(:,:)
1769	DOUBLE PRECISION, ALLOCATABLE :: lon(:,:), lat(:,:), bounds_lon(:,:,:), bounds_lat(:,:,:)
1770	LOGICAL,ALLOCATABLE :: mask(:)
1771	INTEGER :: i,j, ij, n, rank
1772	INTEGER :: nproc_i, nproc_j, nholes, pos_hole
1773	INTEGER,ALLOCATABLE :: size_i(:), begin_i(:), size_j(:), begin_j(:)
1774
1775	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1776	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1777
1778	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
1779	ni_glo=params%ni
1780	nj_glo=params%nj
1781
1782	n=INT(SQRT(mpi_size*1.))
1783	nproc_i=n
1784	nproc_j=n
1785	IF ( n*n == mpi_size) THEN
1786	! do nothing
1787	ELSE IF ( (n+1)*n < mpi_size) THEN
1788	nproc_i=nproc_i+1
1789	nproc_j=nproc_j+1
1790	ELSE
1791	nproc_i=nproc_i+1
1792	ENDIF
1793
1794	nholes=nproc_i*nproc_j-mpi_size
1795
1796	ALLOCATE(size_i(0:nproc_i-1))
1797	ALLOCATE(begin_i(0:nproc_i-1))
1798	DO i=0,nproc_i-1
1799	size_i(i)=ni_glo/nproc_i
1800	IF (i<MOD(ni_glo,nproc_i)) size_i(i)=size_i(i)+1
1801	IF (i==0) THEN
1802	begin_i(i)=0
1803	ELSE
1804	begin_i(i)=begin_i(i-1)+size_i(i-1)
1805	ENDIF
1806	ENDDO
1807
1808	ALLOCATE(size_j(0:nproc_j-1))
1809	ALLOCATE(begin_j(0:nproc_j-1))
1810	DO j=0,nproc_i-1
1811	size_j(j)=nj_glo/nproc_j
1812	IF (j<MOD(nj_glo,nproc_j)) size_j(j)=size_j(j)+1
1813	IF (j==0) THEN
1814	begin_j(j)=0
1815	ELSE
1816	begin_j(j)=begin_j(j-1)+size_j(j-1)
1817	ENDIF
1818	ENDDO
1819
1820
1821	pos_hole=0
1822	rank=0
1823	DO j=0,nproc_j-1
1824	DO i=0,nproc_i-1
1825
1826	ij = j*nproc_i + i
1827	IF (pos_hole<nholes) THEN
1828	IF ( MOD(ij,(nproc_i*nproc_j/nholes)) == 0) THEN
1829	pos_hole=pos_hole+1
1830	CYCLE
1831	ENDIF
1832	ENDIF
1833
1834	IF (mpi_rank==rank) THEN
1835	ibegin = begin_i(i)
1836	ni = size_i(i)
1837	jbegin = begin_j(j)
1838	nj = size_j(j)
1839	ENDIF
1840	rank=rank+1
1841	ENDDO
1842	ENDDO
1843
1844	ALLOCATE(lon_glo(0:ni_glo-1), lat_glo(0:nj_glo-1))
1845	ALLOCATE(bounds_lon_glo(4,0:ni_glo-1), bounds_lat_glo(4,0:nj_glo-1))
1846
1847	DO i=0,ni_glo-1
1848	lon_glo(i)=-180+(i+0.5)*(360./ni_glo)
1849	bounds_lon_glo(1,i)= -180+(i)*(360./ni_glo)
1850	bounds_lon_glo(2,i)= -180+(i+1)*(360./ni_glo)
1851	ENDDO
1852
1853	DO j=0,nj_glo-1
1854	lat_glo(j)=-90+(j+0.5)*(180./nj_glo)
1855	bounds_lat_glo(1,j)= -90+(j)*(180./nj_glo)
1856	bounds_lat_glo(2,j)= -90+(j+1)*(180./nj_glo)
1857	ENDDO
1858
1859	offset_i=2 ! halo of 2 on i
1860	offset_j=1 ! halo of 1 on j
1861
1862	ALLOCATE(lon(0:ni-1,0:nj-1))
1863	ALLOCATE(lat(0:ni-1,0:nj-1))
1864	ALLOCATE(bounds_lon(4,0:ni-1,0:nj-1))
1865	ALLOCATE(bounds_lat(4,0:ni-1,0:nj-1))
1866	! ALLOCATE(mask(0:ni-1,0:nj-1))
1867
1868	ALLOCATE(return_lon(0:ni*nj-1))
1869	ALLOCATE(return_lat(0:ni*nj-1))
1870	ALLOCATE(return_mask(0:ni*nj-1))
1871	ALLOCATE(return_index(0:(ni+2offset_i)(nj+2*offset_j)-1))
1872
1873	ALLOCATE(return_X_lon(0:ni-1))
1874	ALLOCATE(return_X_lat(0:ni-1))
1875	ALLOCATE(return_X_mask(0:ni-1))
1876	ALLOCATE(return_X_index(0:ni-1))
1877
1878	ALLOCATE(return_Y_lon(0:nj-1))
1879	ALLOCATE(return_Y_lat(0:nj-1))
1880	ALLOCATE(return_Y_mask(0:nj-1))
1881	ALLOCATE(return_Y_index(0:nj-1))
1882
1883	return_index=-1
1884	DO j=0,nj-1
1885	DO i=0,ni-1
1886	ij=j*ni+i
1887	return_lon(ij)=lon_glo(ibegin+i)
1888	return_lat(ij)=lat_glo(jbegin+j)
1889	lon(i,j)=return_lon(ij)
1890	lat(i,j)=return_lat(ij)
1891	bounds_lon(1,i,j)=bounds_lon_glo(2,ibegin+i)
1892	bounds_lon(2,i,j)=bounds_lon_glo(1,ibegin+i)
1893	bounds_lon(3,i,j)=bounds_lon_glo(1,ibegin+i)
1894	bounds_lon(4,i,j)=bounds_lon_glo(2,ibegin+i)
1895	bounds_lat(1,i,j)=bounds_lat_glo(1,jbegin+j)
1896	bounds_lat(2,i,j)=bounds_lat_glo(1,jbegin+j)
1897	bounds_lat(3,i,j)=bounds_lat_glo(2,jbegin+j)
1898	bounds_lat(4,i,j)=bounds_lat_glo(2,jbegin+j)
1899
1900	ij=(j+offset_j)(ni+2offset_i)+i+offset_i
1901	return_index(ij)=i+j*ni
1902	ENDDO
1903	ENDDO
1904
1905	CALL set_domain_mask(params, return_lon,return_lat, return_mask)
1906
1907	ALLOCATE(return_X_lon(0:ni-1))
1908	ALLOCATE(return_X_lat(0:ni-1))
1909	ALLOCATE(return_X_mask(0:ni-1))
1910	ALLOCATE(return_X_index(0:ni+2*offset_i-1))
1911
1912	return_X_lat(:)=lat_glo(nj_glo/2)
1913	return_X_lon(:)=lon_glo(ibegin:ibegin+ni-1)
1914	CALL set_domain_mask(params, return_X_lon,return_X_lat, return_X_mask)
1915
1916	return_X_index(:)=-1
1917	DO i=0,ni-1
1918	return_X_index(offset_i+i)=i
1919	ENDDO
1920
1921	ALLOCATE(return_Y_lon(0:nj-1))
1922	ALLOCATE(return_Y_lat(0:nj-1))
1923	ALLOCATE(return_Y_mask(0:nj-1))
1924	ALLOCATE(return_Y_index(0:nj+2*offset_j-1))
1925
1926	return_Y_lat(:)=lat_glo(jbegin:jbegin+nj-1)
1927	return_Y_lon(:)=lon_glo(ni_glo/2)
1928	CALL set_domain_mask(params, return_Y_lon,return_Y_lat, return_Y_mask)
1929
1930	return_Y_index(:)=-1
1931	DO j=0,nj-1
1932	return_Y_index(offset_j+j)=j
1933	ENDDO
1934
1935	return_ni=ni
1936	return_nj=nj
1937
1938	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1939	CALL xios_set_domain_attr(TRIM(domain_id), type="curvilinear", data_dim=2)
1940	CALL xios_set_domain_attr(TRIM(domain_id), ni_glo=ni_glo, ibegin=ibegin, ni=ni, data_ibegin=-offset_i, data_ni=ni+2*offset_i)
1941	CALL xios_set_domain_attr(TRIM(domain_id), nj_glo=nj_glo, jbegin=jbegin, nj=nj, data_jbegin=-offset_j, data_nj=nj+2*offset_j)
1942	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=2, lonvalue_2d=lon, latvalue_2d=lat, mask_1d=return_mask)
1943	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat, nvertex=4)
1944	ENDIF
1945
1946
1947	IF (xios_is_valid_axis(TRIM(domain_id)//"_X")) THEN
1948	CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ni_glo, begin=ibegin, n=ni, data_begin=-offset_i, data_n=ni+2*offset_i, value=return_X_lon, mask=return_X_mask)
1949	! CALL xios_set_axis_attr(TRIM(domain_id)//"_X", n_glo=ni_glo, begin=ibegin, n=ni, data_index=return_X_index, data_n=ni+2*offset_i, value=return_X_lon)
1950	ENDIF
1951
1952	IF (xios_is_valid_axis(TRIM(domain_id)//"_Y")) THEN
1953	CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=nj_glo, begin=jbegin, n=nj, data_begin=-offset_j, data_n=nj+2*offset_j, value=return_Y_lat, mask=return_Y_mask)
1954	! CALL xios_set_axis_attr(TRIM(domain_id)//"_Y", n_glo=nj_glo, begin=jbegin, n=nj, data_index=return_Y_index, data_n=nj+2*offset_j, value=return_Y_lat, mask=return_Y_mask)
1955	ENDIF
1956
1957	END SUBROUTINE init_domain_nemo
1958
1959
1960
1961
1962
1963	SUBROUTINE set_domain_mask(params, lon,lat, mask)
1964	IMPLICIT NONE
1965	TYPE(tmodel_params) :: params
1966	DOUBLE PRECISION :: lon(:)
1967	DOUBLE PRECISION :: lat(:)
1968	LOGICAL :: mask(:)
1969	INTEGER :: i,x
1970
1971	mask(:)=.TRUE.
1972	IF (params%domain_mask) THEN
1973	IF (params%domain_mask_type=="cross") THEN
1974	WHERE (lon(:)-2lat(:)>-10 .AND. lon(:)-2lat(:) <10) mask(:)=.FALSE.
1975	WHERE (2lat(:)+lon(:)>-10 .AND. 2lat(:)+lon(:)<10) mask(:)=.FALSE.
1976	ELSE IF (params%domain_mask_type=="latitude_band") THEN
1977	WHERE (lat(:)>-30 .AND. lat(:)<30) mask(:)=.FALSE.
1978	ENDIF
1979	ENDIF
1980
1981	END SUBROUTINE set_domain_mask
1982
1983
1984	SUBROUTINE set_mask3d(grid_id, params, ni, nj, lon,lat, axis_value)
1985	IMPLICIT NONE
1986	CHARACTER(LEN=*) :: grid_id
1987	TYPE(tmodel_params) :: params
1988	INTEGER :: comm
1989	INTEGER :: ni
1990	INTEGER :: nj
1991	DOUBLE PRECISION, POINTER :: lon(:)
1992	DOUBLE PRECISION, POINTER :: lat(:)
1993	INTEGER, POINTER :: domain_index(:)
1994	DOUBLE PRECISION, POINTER :: axis_value(:)
1995	INTEGER, POINTER :: axis_index(:)
1996	INTEGER ::i,j,ij,k,nk
1997	LOGICAL, ALLOCATABLE :: mask3d(:,:,:)
1998	DOUBLE PRECISION :: r
1999
2000	nk=size(axis_value)
2001
2002	ALLOCATE(mask3d(0:ni-1,0:nj-1,0:nk-1))
2003
2004	mask3d=.TRUE.
2005	DO k=0,nk-1
2006	DO j=0,nj-1
2007	DO i=0,ni-1
2008	ij=j*ni+i
2009	r=sqrt((lon(ij)/2)2 + lat(ij)2) / ((nk-k)*1./nk)
2010	if (r < 60) mask3d(i,j,k)=.FALSE.
2011	ENDDO
2012	ENDDO
2013	ENDDO
2014
2015	IF (params%mask3d) CALL xios_set_grid_attr(grid_id, mask_3d=mask3d)
2016
2017	END SUBROUTINE set_mask3d
2018
2019	SUBROUTINE init_axis_pressure(axis_id,comm,params, return_value, return_mask, return_index)
2020	IMPLICIT NONE
2021	CHARACTER(LEN=*) :: axis_id
2022	TYPE(tmodel_params) :: params
2023	INTEGER :: comm
2024	DOUBLE PRECISION, POINTER :: return_value(:)
2025	LOGICAL, POINTER :: return_mask(:)
2026	INTEGER, POINTER :: return_index(:)
2027
2028	INTEGER :: nlev_glo
2029	INTEGER :: nlev, begin, end
2030	DOUBLE PRECISION, ALLOCATABLE :: value_glo(:)
2031	DOUBLE PRECISION, ALLOCATABLE :: bounds_value_glo(:,:)
2032	DOUBLE PRECISION, ALLOCATABLE :: value(:)
2033	DOUBLE PRECISION, ALLOCATABLE :: bounds_value(:,:)
2034	DOUBLE PRECISION :: dp
2035	INTEGER :: i
2036
2037	INTEGER :: domain_proc_rank, domain_proc_size
2038	INTEGER :: axis_proc_rank, axis_proc_size
2039	INTEGER :: ensemble_proc_size, ensemble_proc_rank
2040
2041	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
2042
2043	nlev_glo=params%nlev
2044
2045	ALLOCATE(value_glo(0:nlev_glo-1), bounds_value_glo(2,0:nlev_glo-1) )
2046
2047	dp=(1-0.1)/(nlev_glo-1)
2048	DO i=0,nlev_glo-1
2049	value_glo(i)=1-i*dp
2050	ENDDO
2051
2052	bounds_value_glo(2,0)=value_glo(0)-(value_glo(1)-value_glo(0))/2
2053	DO i=1,nlev_glo-1
2054	bounds_value_glo(2,i)=(value_glo(i-1)+value_glo(i)) /2
2055	ENDDO
2056
2057	DO i=0,nlev_glo-2
2058	bounds_value_glo(1,i)=(value_glo(i)+value_glo(i+1)) /2
2059	ENDDO
2060	bounds_value_glo(1,nlev_glo-1)=value_glo(nlev_glo-1)-(value_glo(nlev_glo-2)-value_glo(nlev_glo-1))/2
2061
2062	nlev=nlev_glo/axis_proc_size
2063	IF (axis_proc_rank < MOD(nlev_glo,axis_proc_size)) THEN
2064	nlev=nlev+1
2065	begin= axis_proc_rank*nlev
2066	end=begin+nlev-1
2067	ELSE
2068	begin=MOD(nlev_glo,axis_proc_size)(nlev+1) + (axis_proc_rank-MOD(nlev_glo,axis_proc_size))nlev
2069	end=begin+nlev-1
2070	ENDIF
2071
2072	ALLOCATE(value(0:nlev-1), bounds_value(2,0:nlev-1) )
2073	value(:)=value_glo(begin:end)
2074	bounds_value(:,:)=bounds_value_glo(:,begin:end)
2075
2076	ALLOCATE(return_value(0:nlev-1))
2077	ALLOCATE(return_mask(0:nlev-1))
2078	return_value=value
2079	return_mask=.TRUE.
2080	CALL set_axis_mask(params,value,return_mask)
2081	CALL xios_set_axis_attr(axis_id, n_glo=nlev_glo, begin=begin, n=nlev, value=value100000, mask=return_mask, bounds=bounds_value100000, unit='Pa', positive='up')
2082
2083
2084	ALLOCATE(return_index(0:nlev-1))
2085
2086	DO i=0,nlev-1
2087	return_index(i)=i
2088	ENDDO
2089
2090
2091	END SUBROUTINE init_axis_pressure
2092
2093	SUBROUTINE set_axis_mask(params, value, mask)
2094	IMPLICIT NONE
2095	TYPE(tmodel_params) :: params
2096	DOUBLE PRECISION :: value(:)
2097	LOGICAL :: mask(:)
2098	INTEGER :: i,x
2099
2100	mask(:)=.TRUE.
2101	x=size(mask)
2102	IF (params%axis_mask) THEN
2103	DO i=0,x-1
2104	IF (MOD(i,3)==0) mask(i+1)=.FALSE.
2105	IF (MOD(i,4)==0) mask(i+1)=.FALSE.
2106	ENDDO
2107	ENDIF
2108
2109	END SUBROUTINE set_axis_mask
2110
2111	SUBROUTINE init_scalar(scalar_id, comm, params, return_mask)
2112	IMPLICIT NONE
2113	CHARACTER(LEN=*) :: scalar_id
2114	TYPE(tmodel_params) :: params
2115	INTEGER :: comm
2116	LOGICAL :: return_mask
2117	DOUBLE PRECISION :: value =10.
2118
2119	CALL set_scalar_mask(comm, params, return_mask)
2120	CALL xios_set_scalar_attr(scalar_id, value=value, mask=return_mask)
2121
2122	END SUBROUTINE init_scalar
2123
2124	SUBROUTINE set_scalar_mask(comm, params, mask)
2125	IMPLICIT NONE
2126	TYPE(tmodel_params) :: params
2127	INTEGER :: comm
2128	LOGICAL :: mask
2129	INTEGER :: ierr,rank
2130
2131	mask=.TRUE.
2132	IF (params%scalar_mask) THEN
2133	IF (params%scalar_mask_type=="none") THEN
2134	mask=.TRUE.
2135	ELSE IF (params%scalar_mask_type=="full") THEN
2136	mask=.FALSE.
2137	ELSE IF (params%scalar_mask_type=="root") THEN
2138	CALL MPI_COMM_RANK(comm,rank,ierr)
2139	mask = (rank==0)
2140	ELSE IF (params%scalar_mask_type=="sparse") THEN
2141	CALL MPI_COMM_RANK(comm,rank,ierr)
2142	mask = (MOD(rank,2)==0)
2143	ENDIF
2144	ENDIF
2145
2146	END SUBROUTINE set_scalar_mask
2147
2148	SUBROUTINE init_field2D_academic(comm,params, lon, lat, mask, return_field, &
2149	X_lon, X_lat, X_mask, return_fieldX, &
2150	Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
2151	IMPLICIT NONE
2152	TYPE(tmodel_params) :: params
2153	INTEGER :: comm
2154	DOUBLE PRECISION, POINTER :: lon(:)
2155	DOUBLE PRECISION, POINTER :: lat(:)
2156	LOGICAL, POINTER :: mask(:)
2157	DOUBLE PRECISION, POINTER :: return_field(:)
2158
2159	DOUBLE PRECISION, POINTER :: X_lon(:)
2160	DOUBLE PRECISION, POINTER :: X_lat(:)
2161	LOGICAL, POINTER :: X_mask(:)
2162	DOUBLE PRECISION, POINTER :: return_fieldX(:)
2163
2164	DOUBLE PRECISION, POINTER :: Y_lon(:)
2165	DOUBLE PRECISION, POINTER :: Y_lat(:)
2166	LOGICAL, POINTER :: Y_mask(:)
2167	DOUBLE PRECISION, POINTER :: return_fieldY(:)
2168	DOUBLE PRECISION, POINTER :: return_fieldXY(:,:)
2169
2170	DOUBLE PRECISION, PARAMETER :: coef=2., dp_pi=3.14159265359
2171	DOUBLE PRECISION :: dp_length, dp_conv
2172	INTEGER :: i,j,x,y,xy
2173
2174	! Parameters for analytical function
2175	dp_length= 1.2*dp_pi
2176	dp_conv=dp_pi/180.
2177
2178	xy=size(mask)
2179	x=size(X_mask)
2180	y=size(Y_mask)
2181
2182	ALLOCATE(return_field(0:xy-1))
2183
2184	DO i=0,xy-1
2185	IF (mask(i)) THEN
2186	return_field(i)=(coef-SIN(dp_pi(ACOS(COS(lat(i)dp_conv)*&
2187	COS(lon(i)*dp_conv))/dp_length)))
2188	ENDIF
2189	ENDDO
2190
2191
2192	ALLOCATE(return_fieldX(0:x-1))
2193
2194	DO i=0,x-1
2195	IF (X_mask(i)) THEN
2196	return_fieldX(i)=(coef-SIN(dp_pi(ACOS(COS(X_lat(i)dp_conv)*&
2197	COS(X_lon(i)*dp_conv))/dp_length)))
2198	ENDIF
2199	ENDDO
2200
2201
2202	ALLOCATE(return_fieldY(0:y-1))
2203
2204	DO i=0,y-1
2205	IF (Y_mask(i)) THEN
2206	return_fieldY(i)=(coef-SIN(dp_pi(ACOS(COS(Y_lat(i)dp_conv)*&
2207	COS(Y_lon(i)*dp_conv))/dp_length)))
2208	ENDIF
2209	ENDDO
2210
2211	ALLOCATE(return_fieldXY(0:x-1,0:y-1))
2212
2213	DO j=0,y-1
2214	DO i=0,x-1
2215	IF (Y_mask(j) .AND. X_mask(i)) THEN
2216	return_fieldXY(i,j)=(coef-SIN(dp_pi(ACOS(COS(Y_lat(j)dp_conv)*&
2217	COS(X_lon(i)*dp_conv))/dp_length)))
2218	ENDIF
2219	ENDDO
2220	ENDDO
2221
2222	END SUBROUTINE init_field2D_academic
2223
2224
2225	SUBROUTINE init_field2D_constant(comm,params, lon, lat, mask, return_field, &
2226	X_lon, X_lat, X_mask, return_fieldX, &
2227	Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
2228	IMPLICIT NONE
2229	TYPE(tmodel_params) :: params
2230	INTEGER :: comm
2231	DOUBLE PRECISION, POINTER :: lon(:)
2232	DOUBLE PRECISION, POINTER :: lat(:)
2233	LOGICAL, POINTER :: mask(:)
2234	DOUBLE PRECISION, POINTER :: return_field(:)
2235
2236	DOUBLE PRECISION, POINTER :: X_lon(:)
2237	DOUBLE PRECISION, POINTER :: X_lat(:)
2238	LOGICAL, POINTER :: X_mask(:)
2239	DOUBLE PRECISION, POINTER :: return_fieldX(:)
2240
2241	DOUBLE PRECISION, POINTER :: Y_lon(:)
2242	DOUBLE PRECISION, POINTER :: Y_lat(:)
2243	LOGICAL, POINTER :: Y_mask(:)
2244	DOUBLE PRECISION, POINTER :: return_fieldY(:)
2245	DOUBLE PRECISION, POINTER :: return_fieldXY(:,:)
2246	INTEGER :: x,y,xy
2247
2248	xy=size(mask)
2249	x=size(X_mask)
2250	y=size(Y_mask)
2251
2252	ALLOCATE(return_field(0:xy-1))
2253	return_field(:)=1
2254
2255	ALLOCATE(return_fieldX(0:x-1))
2256	return_fieldX=1
2257
2258	ALLOCATE(return_fieldY(0:y-1))
2259	return_fieldY=1
2260
2261	ALLOCATE(return_fieldXY(0:x-1,0:y-1))
2262	return_fieldXY=1
2263
2264	END SUBROUTINE init_field2D_constant
2265
2266	SUBROUTINE init_field2D_rank(comm,params, lon, lat, mask, return_field, &
2267	X_lon, X_lat, X_mask, return_fieldX, &
2268	Y_lon, Y_lat, Y_mask, return_fieldY, return_fieldXY)
2269
2270	IMPLICIT NONE
2271	TYPE(tmodel_params) :: params
2272	INTEGER :: comm
2273	DOUBLE PRECISION, POINTER :: lon(:)
2274	DOUBLE PRECISION, POINTER :: lat(:)
2275	LOGICAL, POINTER :: mask(:)
2276	DOUBLE PRECISION, POINTER :: return_field(:)
2277
2278	DOUBLE PRECISION, POINTER :: X_lon(:)
2279	DOUBLE PRECISION, POINTER :: X_lat(:)
2280	LOGICAL, POINTER :: X_mask(:)
2281	DOUBLE PRECISION, POINTER :: return_fieldX(:)
2282
2283	DOUBLE PRECISION, POINTER :: Y_lon(:)
2284	DOUBLE PRECISION, POINTER :: Y_lat(:)
2285	LOGICAL, POINTER :: Y_mask(:)
2286	DOUBLE PRECISION, POINTER :: return_fieldY(:)
2287	DOUBLE PRECISION, POINTER :: return_fieldXY(:,:)
2288	INTEGER :: rank,ierr
2289	INTEGER :: x,y,xy
2290
2291	CALL MPI_COMM_RANK(comm,rank,ierr)
2292
2293
2294
2295	xy=size(mask)
2296	x=size(X_mask)
2297	y=size(Y_mask)
2298
2299	ALLOCATE(return_field(0:xy-1))
2300	return_field(:)=rank
2301
2302	ALLOCATE(return_fieldX(0:x-1))
2303	return_fieldX=rank
2304
2305	ALLOCATE(return_fieldY(0:y-1))
2306	return_fieldY=rank
2307
2308	ALLOCATE(return_fieldXY(0:x-1,0:y-1))
2309	return_fieldXY=rank
2310
2311	END SUBROUTINE init_field2D_rank
2312
2313
2314
2315	SUBROUTINE get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank, ensemble_proc_size, ensemble_proc_rank)
2316	IMPLICIT NONE
2317	INTEGER,INTENT(IN) :: comm
2318	TYPE(tmodel_params) :: params
2319	INTEGER,INTENT(OUT) :: domain_proc_size
2320	INTEGER,INTENT(OUT) :: domain_proc_rank
2321	INTEGER,INTENT(OUT) :: axis_proc_size
2322	INTEGER,INTENT(OUT) :: axis_proc_rank
2323	INTEGER,INTENT(OUT) :: ensemble_proc_size
2324	INTEGER,INTENT(OUT) :: ensemble_proc_rank
2325
2326	INTEGER :: mpi_rank,mpi_size,rank, ensemble_number
2327	INTEGER :: ierr
2328	INTEGER :: n_domain,n_axis, n_ensemble, min_dist, new_dist, best
2329	INTEGER :: axis_ind, domain_ind
2330	LOGICAL :: axis_layer
2331
2332	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
2333	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
2334
2335	n_ensemble = params%ensemble_proc_n
2336	ensemble_proc_size = mpi_size / n_ensemble
2337	IF ( mpi_rank < MOD(mpi_size,n_ensemble) * (ensemble_proc_size+1) ) THEN
2338	ensemble_proc_size=ensemble_proc_size+1
2339	ensemble_proc_rank = MOD(mpi_rank,ensemble_proc_size)
2340	ensemble_number = mpi_rank / ensemble_proc_size
2341	ELSE
2342	ensemble_number = MOD(mpi_size,n_ensemble)
2343	rank = mpi_rank - MOD(mpi_size,n_ensemble) * (ensemble_proc_size+1)
2344	ensemble_proc_rank= MOD(rank,ensemble_proc_size)
2345	ensemble_number = ensemble_number + rank / ensemble_proc_size
2346	ENDIF
2347
2348	mpi_size=ensemble_proc_size
2349	mpi_rank=ensemble_proc_rank
2350	ensemble_proc_size = n_ensemble
2351	ensemble_proc_rank = ensemble_number
2352
2353	IF (params%axis_proc_n > 0 ) THEN
2354	n_axis=params%axis_proc_n
2355	n_domain = mpi_size / n_axis
2356	axis_layer=.TRUE.
2357	ELSE IF (params%domain_proc_n > 0 ) THEN
2358	n_domain=params%domain_proc_n
2359	n_axis = mpi_size / n_domain
2360	axis_layer=.FALSE.
2361	ELSE
2362	IF (params%axis_proc_frac==0) THEN
2363	params%axis_proc_frac=1
2364	params%domain_proc_frac=mpi_size
2365	ELSE IF (params%domain_proc_frac==0) THEN
2366	params%domain_proc_frac=1
2367	params%axis_proc_frac=mpi_size
2368	ENDIF
2369
2370	n_domain = INT(sqrt(params%domain_proc_frac * mpi_size/ params%axis_proc_frac))
2371	n_axis = INT(sqrt(params%axis_proc_frac * mpi_size/ params%domain_proc_frac))
2372
2373
2374	min_dist= mpi_size - n_domain*n_axis
2375	best=0
2376
2377	new_dist = mpi_size -(n_domain+1)*n_axis
2378	IF (new_dist < min_dist .AND. new_dist >= 0 ) THEN
2379	min_dist=new_dist
2380	best=1
2381	ENDIF
2382
2383	new_dist=mpi_size-n_domain*(n_axis+1)
2384	IF (new_dist < min_dist .AND. new_dist >= 0 ) THEN
2385	min_dist=new_dist
2386	best=2
2387	ENDIF
2388
2389	IF (best==0) THEN
2390	ELSE IF (best==1) THEN
2391	n_domain=n_domain+1
2392	ELSE IF (best==2) THEN
2393	n_axis=n_axis+1
2394	ENDIF
2395
2396	IF ( MOD(mpi_size,n_axis) <= MOD(mpi_size,n_domain)) axis_layer=.TRUE.
2397
2398	ENDIF
2399
2400	IF ( axis_layer) THEN
2401	!!! n_axis layer
2402	IF (mpi_rank < MOD(mpi_size,n_axis)*(n_domain+1)) THEN
2403	axis_proc_rank=mpi_rank/(n_domain+1)
2404	domain_proc_rank=MOD(mpi_rank,n_domain+1)
2405	axis_proc_size=n_axis
2406	domain_proc_size=n_domain+1
2407	ELSE
2408	rank=mpi_rank-MOD(mpi_size,n_axis)*(n_domain+1)
2409	axis_proc_size=n_axis
2410	axis_proc_rank=MOD(mpi_size,n_axis)+rank/n_domain
2411	domain_proc_rank=MOD(rank,n_domain)
2412	domain_proc_size=n_domain
2413	ENDIF
2414	ELSE
2415	!!! n_domain column
2416	IF (mpi_rank < MOD(mpi_size,n_domain)*(n_axis+1)) THEN
2417	domain_proc_rank=mpi_rank/(n_axis+1)
2418	axis_proc_rank=MOD(mpi_rank,n_axis+1)
2419	domain_proc_size=n_domain
2420	axis_proc_size=n_axis+1
2421	ELSE
2422	rank=mpi_rank-MOD(mpi_size,n_domain)*(n_axis+1)
2423	domain_proc_size=n_domain
2424	domain_proc_rank=MOD(mpi_size,n_domain)+rank/n_axis
2425	axis_proc_rank=MOD(rank,n_axis)
2426	axis_proc_size=n_axis
2427	ENDIF
2428	ENDIF
2429
2430
2431	END SUBROUTINE get_decomposition
2432
2433	END PROGRAM generic_testcase
2434
2435

Note: See TracBrowser for help on using the repository browser.

Download in other formats: