source: XIOS/dev/dev_trunk_graph/src/test/generic_testcase.f90 @ 2027

Last change on this file since 2027 was 2027, checked in by yushan, 9 months ago

Graph intermediate commit to a tmp branch. Integrate latest modifications of branch coupling

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