source: XIOS/trunk/src/test/generic_testcase.f90 @ 1921

Last change on this file since 1921 was 1921, checked in by yushan, 2 years ago

trunk : complete axis transformation tests

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