Context Navigation

source: XIOS3/trunk/src/test/generic_testcase.f90 @ 2532

Last change on this file since 2532 was 2520, checked in by jderouillat, 13 months ago
Replace MPI probing on intercommunicator by probing on intracommunicator
Property svn:executable set to ``*
File size: 89.8 KB

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

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

Download in other formats: