Context Navigation

source: XIOS/trunk/src/test/generic_testcase.f90 @ 1921

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

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

Download in other formats: