Context Navigation

source: XIOS/dev/dev_trunk_omp/src/test/generic_testcase.f90 @ 1838

Last change on this file since 1838 was 1726, checked in by yushan, 5 years ago
Generic_testcase : Cleanup and better organization
File size: 84.0 KB

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

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

Download in other formats: