source: XIOS/xios_training/hands-on-11/TP_src/generic_testcase.f90 @ 2102

Last change on this file since 2102 was 2102, checked in by yushan, 15 months ago

XIOS TRAINING add Tp source for HO 11

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