Context Navigation

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

Last change on this file since 1674 was 1674, checked in by ymipsl, 5 years ago

Add generic testcase which would replace all previous testcase.

YM

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

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

Download in other formats: