Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

generic_testcase.f90

Last change on this file was 2591, checked in by jderouillat, 8 months ago
Free additional communicators
Property svn:executable set to ``*
File size: 89.9 KB

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

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

Context Navigation

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

Download in other formats: