Context Navigation

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

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

Fix problem in remoteConnector when computing grid to sent to server.
Some optimisations when grid is not distributed need knowledge of the workflow view.
New CGridClientServerConnector class created based on CGridRemoteConnector.

YM

Property svn:executable set to *

File size: 87.2 KB

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

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

Download in other formats: