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

Last change on this file since 1926 was 1926, checked in by ymipsl, 14 months ago

Bug fix in generic testcase for orchidee grid.

YM

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