Context Navigation

source: XIOS2/trunk/src/test/generic_testcase.f90 @ 2442

Last change on this file since 2442 was 2442, checked in by jderouillat, 19 months ago
Fix bounds management in transformations. Disabled temporarily bounds settings in the rectilinear case of the generic_testcase
File size: 86.3 KB

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

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

Download in other formats: