Context Navigation

source: XIOS2/trunk/src/test/generic_testcase.f90 @ 2440

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

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

Download in other formats: