Context Navigation

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

Last change on this file since 2280 was 1980, checked in by yushan, 4 years ago
trunk : axis interpolate can have coordinate source (coordinate_src) and coordinate destination (coordinate_dst), while previous attribute coordinate compatible to source
File size: 85.2 KB

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

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

Download in other formats: