Context Navigation

source: XIOS/trunk/src/test/test_interpolate.f90 @ 1980

Last change on this file since 1980 was 1980, checked in by yushan, 4 years ago
trunk : axis interpolate can have coordinate source (coordinate_src) and coordinate destination (coordinate_dst), while previous attribute coordinate compatible to source
File size: 85.5 KB

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

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

Download in other formats: