Context Navigation

source: XIOS/dev/dev_trunk_omp/src/test/generic_testcase.f90 @ 1702

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

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

Download in other formats: