Context Navigation

source: XIOS/dev/dev_ym/XIOS_COUPLING/src/test/generic_testcase.f90 @ 2175

Last change on this file since 2175 was 2175, checked in by ymipsl, 3 years ago

For generic testcase, for lmdz domain type, mask was not set for X and Y derived axis.

YM

Property svn:executable set to *

File size: 87.2 KB

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

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

Download in other formats: