Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

oasis_testcase.f90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 15 months ago
Moved oasis-mct_5.0 in oasis3-mct/branches directory.
File size: 57.5 KB

Line
1	PROGRAM oasis_testcase
2	USE xios
3	USE mod_wait
4	USE netcdf
5	IMPLICIT NONE
6	INCLUDE "mpif.h"
7	INTEGER,PARAMETER :: unit=10
8	INTEGER,PARAMETER :: len_str=255
9	INTEGER :: ierr
10
11	INTEGER :: nb_proc_toy
12	CHARACTER(LEN=len_str) :: start_date
13	CHARACTER(LEN=len_str) :: duration
14	NAMELIST /params_run/ nb_proc_toy, duration
15
16
17	TYPE tmodel_params
18	CHARACTER(len_str) :: grids_file=""
19	CHARACTER(len_str) :: masks_file=""
20	CHARACTER(len_str) :: areas_file=""
21	CHARACTER(len_str) :: timestep=""
22	CHARACTER(len_str) :: domain=""
23	LOGICAL :: is_ocean
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 :: ni
30	INTEGER :: nj
31	CHARACTER(len_str) :: init_field2D=""
32	DOUBLE PRECISION :: gf_coef
33	LOGICAL :: domain_mask
34	END TYPE tmodel_params
35
36	LOGICAL :: i_am_atm, i_am_server
37	INTEGER :: rank, size_loc
38
39
40	OPEN(unit=unit, file='param.def',status='old',iostat=ierr)
41	nb_proc_toy=1
42	duration="1h"
43	start_date="2018-01-01"
44	READ (unit, nml=params_run)
45	CLOSE(unit)
46
47	! MPI Initialization
48
49	CALL MPI_INIT(ierr)
50	CALL init_wait
51	CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
52	CALL MPI_COMM_SIZE(MPI_COMM_WORLD,size_loc,ierr)
53	IF (size_loc < nb_proc_toy) THEN
54	PRINT *,"inconsistency : size=",size_loc," nb_proc_toy=",nb_proc_toy
55	CALL MPI_ABORT()
56	ENDIF
57
58	i_am_atm=.FALSE. ; i_am_server=.FALSE.
59
60	! First ranks are dedicated to Toy last to Server
61
62	IF (rank < nb_proc_toy) THEN
63	i_am_atm=.TRUE.
64	ELSE
65	i_am_server=.TRUE.
66	ENDIF;
67
68
69	IF (i_am_server) THEN
70	CALL xios_init_server()
71	ELSE
72	IF (i_am_atm) CALL model("toy")
73	CALL xios_finalize()
74	ENDIF
75
76	WRITE(*,'(A,I0,A)') "Rank ",rank, " finished Successfully"
77	CALL MPI_FINALIZE(ierr)
78
79	CONTAINS
80
81	SUBROUTINE model(model_id)
82	IMPLICIT NONE
83	CHARACTER(LEN=*) :: model_id
84	TYPE(tmodel_params) :: params, other_params
85	TYPE(xios_context) :: ctx_hdl
86	INTEGER :: comm
87	TYPE(xios_date) :: cdate, edate
88	TYPE(xios_duration) :: dtime
89	INTEGER :: rank
90	INTEGER :: size_loc
91	INTEGER :: ts
92
93	DOUBLE PRECISION, POINTER :: lon(:)
94	DOUBLE PRECISION, POINTER :: lat(:)
95	LOGICAL, POINTER :: domain_mask(:)
96	INTEGER, POINTER :: domain_index(:)
97
98	DOUBLE PRECISION, POINTER :: field2D_init(:)
99
100	DOUBLE PRECISION, POINTER :: field2D(:), other_field2D(:)
101
102	INTEGER :: ni,nj
103	INTEGER :: i,xy
104	INTEGER :: ierr
105
106	LOGICAL :: ok_field2D
107	LOGICAL :: ok_other_field2D
108
109	!! XIOS Initialization (get the local communicator)
110	CALL xios_initialize(TRIM(model_id),return_comm=comm)
111	CALL MPI_COMM_RANK(comm,rank,ierr)
112	CALL MPI_COMM_SIZE(comm,size_loc,ierr)
113	CALL xios_context_initialize(TRIM(model_id),comm)
114	CALL xios_get_handle(TRIM(model_id),ctx_hdl)
115	CALL xios_set_current_context(ctx_hdl)
116
117	CALL init_model_params('src_',params)
118	CALL init_model_params('tgt_',other_params)
119
120	IF (params%is_ocean .OR. other_params%is_ocean) THEN
121	params%domain_mask = .TRUE.
122	other_params%domain_mask = .TRUE.
123	END IF
124
125	IF (TRIM(params%init_field2D) == "genmask") THEN
126	other_params%domain_mask = .FALSE.
127	IF (.NOT. params%is_ocean) THEN
128	IF (rank == 0) WRITE(,) &
129	& '====> WARNING: meaningless fractionary mask computation if the source mesh is not oceanic'
130	END IF
131	END IF
132
133	!!! Definition de la start date et du timestep
134	CALL xios_set_start_date(xios_date_convert_from_string(start_date//" 00:00:00"))
135	dtime=xios_duration_convert_from_string(TRIM(params%timestep))
136	CALL xios_set_timestep(timestep=dtime)
137
138	!!! Calcul de temps elaps par seconde pour respecter le SYPD (hyp : pas de delai d'I/O)
139	CALL xios_get_start_date(cdate)
140	edate=cdate+xios_duration_convert_from_string(TRIM(duration))
141
142	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143	!!! standard initialisation of domain, grid, field !
144	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145
146	CALL init_domain("src_domain", comm, params, ni, nj, lon, lat, domain_mask, domain_index)
147
148	CALL init_field2D(comm, params, lon, lat, domain_mask, field2D_init)
149
150	xy=SIZE(domain_index)
151
152	ALLOCATE(field2D(0:xy-1))
153
154	DO i=0,xy-1
155	IF (domain_index(i)/=-1) THEN
156	field2D(i)=field2D_init(domain_index(i))
157	ENDIF
158	ENDDO
159
160	ok_field2D=xios_is_valid_field("src_field2D") ;
161
162	IF (ASSOCIATED(lon)) DEALLOCATE(lon)
163	IF (ASSOCIATED(lat)) DEALLOCATE(lat)
164	IF (ASSOCIATED(domain_mask)) DEALLOCATE(domain_mask)
165	IF (ASSOCIATED(domain_index)) DEALLOCATE(domain_index)
166	IF (ASSOCIATED(field2D_init)) DEALLOCATE(field2D_init)
167
168	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
169	! duplicated section for other_ !
170	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171
172	CALL init_domain("tgt_domain", comm, other_params, ni, nj, lon, lat, domain_mask, domain_index)
173
174	CALL init_field2D(comm, other_params, lon, lat, domain_mask, field2D_init)
175
176	xy=SIZE(domain_index)
177
178	ALLOCATE(other_field2D(0:xy-1))
179
180	DO i=0,xy-1
181	IF (domain_index(i)/=-1) THEN
182	other_field2D(i)=field2D_init(domain_index(i))
183	ENDIF
184	ENDDO
185
186	ok_other_field2D=xios_is_valid_field("tgt_field2D") ;
187
188
189	!!!!!!!!!!!!!!!!!!!!! end of duplicated section !!!!!!!!!!!!!!!!!!!!!!!!!!
190
191	! Conditional setting of the interpolation for the fractionary mask
192
193	IF (TRIM(params%init_field2D) == "genmask") THEN
194	CALL xios_set_interpolate_domain_attr("interpolation", &
195	& order=1, renormalize=.FALSE., use_area=.FALSE., &
196	& detect_missing_value=.FALSE., write_weight=.FALSE.)
197	END IF
198
199	!
200
201	ts=1
202	cdate=cdate+dtime
203	CALL xios_close_context_definition()
204	CALL xios_set_current_context(TRIM(model_id))
205
206
207	DO WHILE ( cdate <= edate )
208
209	!!! Mise a jour du pas de temps
210	CALL xios_update_calendar(ts)
211
212	IF (ok_field2D) CALL xios_send_field("src_field2D",field2D)
213
214	field2D=field2D+1
215
216	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
217	! duplicated section for other_ !
218	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
219
220
221	IF (ok_other_field2D) CALL xios_send_field("tgt_field2D", other_field2D)
222	other_field2D=other_field2D+1
223
224
225	!!!!!!!!!!!!!!!!!!!!!! end of duplicated section !!!!!!!!!!!!!!!!!
226
227
228	cdate=cdate+dtime
229	ts=ts+1
230	ENDDO
231
232	CALL xios_context_finalize()
233	CALL MPI_COMM_FREE(comm, ierr)
234
235	END SUBROUTINE model
236
237	SUBROUTINE init_model_params(prefix,params)
238	IMPLICIT NONE
239	CHARACTER(LEN=*) :: prefix
240	TYPE(tmodel_params) :: params
241
242
243	IF (.NOT. xios_getvar(prefix//"timestep", params%timestep) ) params%timestep="1h"
244	IF (.NOT. xios_getvar(prefix//"domain", params%domain) ) params%domain="lmdz"
245	IF (.NOT. xios_getvar("grids_input_file", params%grids_file) ) params%grids_file="input_data/grids.nc"
246	IF (.NOT. xios_getvar("masks_input_file", params%masks_file) ) params%masks_file="input_data/masks.nc"
247	IF (.NOT. xios_getvar("areas_input_file", params%areas_file) ) params%areas_file="input_data/areas.nc"
248	IF (.NOT. xios_getvar(prefix//"domain_mask", params%domain_mask) ) params%domain_mask=.FALSE.
249	IF (.NOT. xios_getvar(prefix//"ni", params%ni) ) params%ni=36
250	IF (.NOT. xios_getvar(prefix//"nj", params%nj) ) params%nj=18
251	IF (.NOT. xios_getvar("init_field2D", params%init_field2D) ) params%init_field2D="academic"
252	IF (.NOT. xios_getvar("gulf_stream_coeff", params%gf_coef) ) params%gf_coef=0.0
253
254	IF (.NOT. xios_getvar(prefix//"domain_proc_n", params%domain_proc_n)) params%domain_proc_n=0
255	IF (.NOT. xios_getvar(prefix//"axis_proc_n", params%axis_proc_n)) params%axis_proc_n=0
256	IF ((.NOT. xios_getvar(prefix//"domain_proc_frac", params%domain_proc_frac)) .AND. &
257	(.NOT. xios_getvar(prefix//"axis_proc_frac", params%axis_proc_frac))) THEN
258	params%domain_proc_frac=1.0
259	params%axis_proc_frac=0.0
260	ELSE IF (.NOT. xios_getvar(prefix//"domain_proc_frac", params%domain_proc_frac)) THEN
261	params%domain_proc_frac=0.0
262	ELSE IF (.NOT. xios_getvar(prefix//"axis_proc_frac", params%axis_proc_frac)) THEN
263	params%axis_proc_frac=0.0
264	ENDIF
265
266	SELECT CASE(TRIM(params%domain))
267	CASE('torc', 'nogt', 'to25', 't12e' )
268	params%is_ocean = .TRUE.
269	CASE DEFAULT
270	params%is_ocean = .FALSE.
271	END SELECT
272
273	params%domain_mask = params%domain_mask .OR. params%is_ocean
274
275	END SUBROUTINE init_model_params
276
277
278	SUBROUTINE init_domain(domain_id, comm, params, return_ni, return_nj, &
279	return_lon, return_lat, return_mask, return_index)
280	IMPLICIT NONE
281	CHARACTER(LEN=*),INTENT(IN) :: domain_id
282	TYPE(tmodel_params),INTENT(IN) :: params
283	INTEGER,INTENT(IN) :: comm
284	INTEGER,INTENT(OUT) :: return_ni
285	INTEGER,INTENT(OUT) :: return_nj
286	DOUBLE PRECISION, POINTER :: return_lon(:)
287	DOUBLE PRECISION, POINTER :: return_lat(:)
288	LOGICAL, POINTER :: return_mask(:)
289	INTEGER, POINTER :: return_index(:)
290
291	SELECT CASE(params%domain)
292	CASE("lmdz")
293	CALL init_domain_lmdz(domain_id,comm,params, return_ni, return_nj, &
294	return_lon, return_lat, return_mask, return_index)
295	CASE("bggd", "tlan")
296	CALL init_domain_bggd(domain_id,comm,params, return_ni, return_nj, &
297	return_lon, return_lat, return_mask, return_index)
298	CASE("icos", "icoh")
299	CALL init_domain_icos(domain_id,comm,params, return_ni, return_nj, &
300	return_lon, return_lat, return_mask, return_index)
301	CASE("ssea", "t359", "sse7")
302	CALL init_domain_ssea(domain_id,comm,params, return_ni, return_nj, &
303	return_lon, return_lat, return_mask, return_index)
304	CASE("gaussian")
305	CALL init_domain_gaussian(domain_id,comm,params, return_ni, return_nj, &
306	return_lon, return_lat, return_mask, return_index)
307	CASE("nemo")
308	CALL init_domain_nemo(domain_id,comm,params, return_ni, return_nj, &
309	return_lon, return_lat, return_mask, return_index)
310	CASE("nogt", "torc", "to25", "t12e")
311	CALL init_domain_orca(domain_id,comm,params, return_ni, return_nj, &
312	return_lon, return_lat, return_mask, return_index)
313	END SELECT
314
315	END SUBROUTINE init_domain
316
317	SUBROUTINE init_field2D(comm, params, lon, lat, mask, return_field)
318	IMPLICIT NONE
319	INTEGER :: comm
320	TYPE(tmodel_params) :: params
321	DOUBLE PRECISION, POINTER :: lon(:)
322	DOUBLE PRECISION, POINTER :: lat(:)
323	LOGICAL, POINTER :: mask(:)
324	DOUBLE PRECISION, POINTER :: return_field(:)
325
326	SELECT CASE(params%init_field2D)
327	CASE("academic")
328	CALL init_field2D_academic(params,lon, lat, mask, return_field)
329	CASE("vortex")
330	CALL init_field2D_vortex(params,lon, lat, mask, return_field)
331	CASE("one")
332	CALL init_field2D_one(mask, return_field)
333	CASE("rank")
334	CALL init_field2D_rank(comm, mask, return_field)
335	CASE("genmask")
336	CALL init_field2D_mask(comm, mask, return_field)
337	END SELECT
338
339	END SUBROUTINE init_field2D
340
341
342	SUBROUTINE init_domain_gaussian(domain_id, comm, params, return_ni, return_nj, &
343	return_lon,return_lat,return_mask, return_index)
344	IMPLICIT NONE
345	CHARACTER(LEN=*) :: domain_id
346	TYPE(tmodel_params) :: params
347	INTEGER :: comm
348	INTEGER :: return_ni
349	INTEGER :: return_nj
350	DOUBLE PRECISION, POINTER :: return_lon(:)
351	DOUBLE PRECISION, POINTER :: return_lat(:)
352	LOGICAL, POINTER :: return_mask(:)
353	INTEGER, POINTER :: return_index(:)
354	INTEGER :: domain_proc_rank, domain_proc_size
355
356	INTEGER :: mpi_rank, mpi_size
357	INTEGER :: ierr
358
359	INTEGER :: nlon, nlat
360	DOUBLE PRECISION,ALLOCATABLE :: lon_glo(:),lat_glo(:)
361	DOUBLE PRECISION,ALLOCATABLE :: bounds_lon_glo(:,:),bounds_lat_glo(:,:)
362	INTEGER,ALLOCATABLE :: i_index_glo(:)
363	INTEGER,ALLOCATABLE :: i_index(:)
364	LOGICAL,ALLOCATABLE :: mask_glo(:),mask(:)
365	DOUBLE PRECISION,ALLOCATABLE :: lon(:),lat(:),field_A_srf(:,:)
366	DOUBLE PRECISION,ALLOCATABLE :: bounds_lon(:,:),bounds_lat(:,:)
367
368	INTEGER :: i,j,n
369	INTEGER :: ncell_glo,ncell,ind
370	REAL :: ilon,ilat
371	DOUBLE PRECISION, PARAMETER :: Pi=3.14159265359
372	INTEGER,ALLOCATABLE :: list_ind(:,:)
373	INTEGER :: rank,j1,j2,np,ncell_x
374	INTEGER :: data_n_index
375	INTEGER,ALLOCATABLE :: data_i_index(:)
376	INTEGER :: axis_proc_rank, axis_proc_size
377
378
379
380	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
381	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
382
383	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank)
384
385	mpi_rank=domain_proc_rank
386	mpi_size=domain_proc_size
387	nlon=params%ni
388	nlat=params%nj
389
390
391	ncell_glo=0
392	DO j=1,nlat
393	n = NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
394	IF (n<8) n=8
395	ncell_glo=ncell_glo+n
396	ENDDO
397
398	ALLOCATE(lon_glo(ncell_glo))
399	ALLOCATE(lat_glo(ncell_glo))
400	ALLOCATE(bounds_lon_glo(4,ncell_glo))
401	ALLOCATE(bounds_lat_glo(4,ncell_glo))
402	ALLOCATE(i_index_glo(ncell_glo)) ; i_index_glo(:) = -1
403	ALLOCATE(mask_glo(ncell_glo))
404	ALLOCATE(list_ind(nlon,nlat))
405
406	ind=0
407	DO j=1,nlat
408	n = NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
409	IF (mpi_rank == 0 .AND. j==1) PRINT*,"--- Gauss lon at first lat ",n
410	IF (mpi_rank == 0 .AND. j==nlat) PRINT*,"--- Gauss lon at last lat ",n
411	IF (n<8) n=8
412
413	DO i=1,n
414	ind=ind+1
415	list_ind(i,j)=ind
416	ilon=i-0.5
417	ilat=j-0.5
418
419	lat_glo(ind)= 90-(ilat*180./nlat)
420	lon_glo(ind)= (ilon*360./n)-180.
421
422
423	bounds_lat_glo(1,ind)= 90-((ilat-0.5)*180./nlat)
424	bounds_lon_glo(1,ind)=((ilon-0.5)*360./n)-180.
425
426	bounds_lat_glo(2,ind)= 90-((ilat-0.5)*180./nlat)
427	bounds_lon_glo(2,ind)=((ilon+0.5)*360./n)-180.
428
429	bounds_lat_glo(3,ind)= 90-((ilat+0.5)*180./nlat)
430	bounds_lon_glo(3,ind)=((ilon+0.5)*360./n)-180.
431
432	bounds_lat_glo(4,ind)= 90-((ilat+0.5)*180./nlat)
433	bounds_lon_glo(4,ind)=((ilon-0.5)*360./n)-180.
434
435	ENDDO
436	ENDDO
437
438	! mpi_size=32
439	rank=(mpi_size-1)/2
440	ncell_x=SQRT(ncell_glo*1./mpi_size)
441
442	j1=nlat/2
443	DO WHILE(rank>=0)
444	j2=MAX(j1-ncell_x+1,1)
445	j=(j1+j2)/2
446	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
447	IF (n<8) n=8
448	np = MIN(n/ncell_x,rank+1) ;
449	IF (j2==1) np=rank+1
450	IF (np == rank+1) j2 = 1
451
452	IF (mpi_rank == 0 ) PRINT *,"domain ",j2,j1,rank,np ;
453	DO j=j2,j1
454	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
455	IF (n<8) n=8
456	DO i=1,n
457	ind=list_ind(i,j)
458	IF ( (i-1) < MOD(n,np)*(n/np+1)) THEN
459	i_index_glo(ind) = rank - (i-1)/(n/np+1)
460	ELSE
461	i_index_glo(ind) = rank-(MOD(n,np)+ (i-1-MOD(n,np)*(n/np+1))/(n/np))
462	ENDIF
463	ENDDO
464	ENDDO
465	rank=rank-np
466	j1=j2-1
467	ENDDO
468
469	rank=MIN((mpi_size-1)/2+1, (mpi_size-1))
470	ncell_x=SQRT(ncell_glo*1./mpi_size)
471
472	j1=nlat/2+1
473	DO WHILE(rank<=mpi_size-1)
474	j2=MIN(j1+ncell_x-1,nlat)
475	j=(j1+j2)/2
476	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
477	IF (n<8) n=8
478	np = MIN(n/ncell_x,mpi_size-rank) ;
479	IF (j2==nlat) np=mpi_size-rank
480	IF (np==mpi_size-rank) j2 = nlat
481
482	IF (mpi_rank == 0 ) PRINT *,"domain ",j2,j1,rank,np ;
483	DO j=j1,j2
484	n=NINT(COS(Pi/2-(j-0.5)PI/nlat)nlon)
485	IF (n<8) n=8
486	DO i=1,n
487	ind=list_ind(i,j)
488	IF ( (i-1) < MOD(n,np)*(n/np+1)) THEN
489	i_index_glo(ind) = rank + (i-1)/(n/np+1)
490	ELSE
491	i_index_glo(ind) = rank+(MOD(n,np)+ (i-1-MOD(n,np)*(n/np+1))/(n/np))
492	ENDIF
493	ENDDO
494	ENDDO
495	rank=rank+np
496	j1=j2+1
497	ENDDO
498
499	ncell=0
500	DO ind=1,ncell_glo
501	IF (i_index_glo(ind)==mpi_rank) ncell=ncell+1
502	ENDDO
503	ALLOCATE(i_index(ncell))
504	ALLOCATE(lon(ncell))
505	ALLOCATE(lat(ncell))
506	ALLOCATE(bounds_lon(4,ncell))
507	ALLOCATE(bounds_lat(4,ncell))
508	ALLOCATE(mask(ncell))
509	ncell=0
510	data_n_index=0
511	DO ind=1,ncell_glo
512	IF (i_index_glo(ind)==mpi_rank) THEN
513	ncell=ncell+1
514	i_index(ncell)=ind-1
515	lon(ncell)=lon_glo(ind)
516	lat(ncell)=lat_glo(ind)
517	bounds_lon(:,ncell)=bounds_lon_glo(:,ind)
518	bounds_lat(:,ncell)=bounds_lat_glo(:,ind)
519	field_A_srf(ncell,:)=i_index_glo(ind)
520	IF (MOD(ind,8)>=0 .AND. MOD(ind,8)<2) THEN
521	mask(ncell)=.TRUE.
522	data_n_index=data_n_index+1
523	ELSE
524	mask(ncell)=.TRUE.
525	data_n_index=data_n_index+1
526	ENDIF
527	ENDIF
528	ENDDO
529
530	ALLOCATE(data_i_index(data_n_index))
531	data_n_index=0
532	DO ind=1,ncell
533	IF (mask(ind)) THEN
534	data_n_index=data_n_index+1
535	data_i_index(data_n_index)=ind-1
536	ENDIF
537	ENDDO
538
539	ALLOCATE(return_lon(0:ncell-1))
540	ALLOCATE(return_lat(0:ncell-1))
541	ALLOCATE(return_mask(0:ncell-1))
542	ALLOCATE(return_index(0:data_n_index-1))
543	return_lon(0:ncell-1)=lon(1:ncell)
544	return_lat(0:ncell-1)=lat(1:ncell)
545	return_index(0:data_n_index-1)= data_i_index(1:data_n_index)
546	CALL set_domain_mask(params, lon, lat, return_mask)
547
548	return_ni=ncell
549	return_nj=1
550
551
552	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
553	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="unstructured", ni_glo=ncell_glo, ni=ncell, ibegin=0, i_index=i_index)
554	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)
555	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_1D=lon, latvalue_1D=lat, nvertex=4, bounds_lon_1D=bounds_lon, &
556	bounds_lat_1D=bounds_lat)
557	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
558	ENDIF
559
560	END SUBROUTINE init_domain_gaussian
561
562	SUBROUTINE init_domain_ssea(domain_id, comm, params, return_ni, return_nj, &
563	return_lon,return_lat,return_mask, return_index)
564	IMPLICIT NONE
565	CHARACTER(LEN=*) :: domain_id
566	TYPE(tmodel_params) :: params
567
568	INTEGER :: comm
569	INTEGER :: return_ni
570	INTEGER :: return_nj
571	DOUBLE PRECISION, POINTER :: return_lon(:)
572	DOUBLE PRECISION, POINTER :: return_lat(:)
573	LOGICAL, POINTER :: return_mask(:)
574	INTEGER, POINTER :: return_index(:)
575	INTEGER :: mpi_rank, mpi_size
576	INTEGER :: ierr
577
578	DOUBLE PRECISION, POINTER :: lon_glo2d(:,:), lat_glo2d(:,:)
579	DOUBLE PRECISION, POINTER :: clo_glo2d(:,:,:), cla_glo2d(:,:,:)
580	INTEGER, POINTER :: msk_glo2d(:,:)
581	DOUBLE PRECISION, POINTER :: srf_glo2d(:,:)
582
583	INTEGER :: ni_glo,nj_glo, ncrn, nbp
584
585	INTEGER :: i, j, ipx, ic, nlat, npx, npy, ib, ie
586	INTEGER, ALLOCATABLE :: glo_rank(:)
587	INTEGER, ALLOCATABLE :: ila_pcx(:)
588	DOUBLE PRECISION :: last_lat, curr_lon
589	INTEGER, ALLOCATABLE :: i_index(:)
590	DOUBLE PRECISION, ALLOCATABLE :: srf_loc2d(:,:)
591	DOUBLE PRECISION, ALLOCATABLE :: bounds_lon(:,:), bounds_lat(:,:)
592
593	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
594	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
595
596	CALL read_oasis_grid(params, ni_glo, nj_glo, ncrn,&
597	& lon_glo2d, lat_glo2d, clo_glo2d, cla_glo2d, msk_glo2d, srf_glo2d)
598
599	! Count the latitudes
600	nlat = 1
601	last_lat = lat_glo2d(0,0)
602	DO i = 1, ni_glo-1
603	IF (ABS(lat_glo2d(i,0) - last_lat) > 2.*EPSILON(1.0)) THEN
604	nlat = nlat + 1
605	last_lat = lat_glo2d(i,0)
606	END IF
607	END DO
608
609	! Choose the partition
610	npy=INT(SQRT(mpi_size*1.))
611	DO WHILE(MOD(mpi_size,npy) /= 0)
612	npy = npy - 1
613	END DO
614	npx = mpi_size / npy
615
616	! Decomposition
617	ALLOCATE(glo_rank(0:ni_glo-1))
618	ALLOCATE(ila_pcx(0:npx-1))
619	ib = 0
620	DO i = 0, npy-1
621	! First step: decompose the domain in pseudo latitude bands
622	ie = ib + ni_glo/npy
623	IF (i < MOD(ni_glo,npy)) ie = ie + 1
624	glo_rank(ib:ie-1) = i*npx
625	! Second step: decompose the pseudo latitude bands
626	ila_pcx(:) = 0
627	DO j = ib, ie-1
628	curr_lon = lon_glo2d(j,0)
629	IF (curr_lon > 360.0) curr_lon = curr_lon - 360.
630	IF (curr_lon < 0.0) curr_lon = curr_lon + 360.
631	ipx = MIN(FLOOR(npx*curr_lon/360.),npx-1)
632	glo_rank(j) = glo_rank(j) + ipx
633	ila_pcx(ipx) = ila_pcx(ipx) + 1
634	END DO
635	! Third step: balance them
636	DO ipx = 0, npx - 1
637	ila_pcx(ipx) = ila_pcx(ipx) - (ie-ib)/npx
638	IF (ipx < MOD(ie-ib,npx)) ila_pcx(ipx) = ila_pcx(ipx) - 1
639	END DO
640	DO ipx = 0, npx - 2
641	DO WHILE (ila_pcx(ipx) /=0)
642	j = ib
643	DO WHILE (j < ie-1)
644	IF (glo_rank(j) == i*npx + ipx .AND. &
645	& glo_rank(j) == glo_rank(j+1)-1) THEN
646	IF (ila_pcx(ipx) > 0) THEN
647	glo_rank(j) = glo_rank(j+1)
648	ila_pcx(ipx) = ila_pcx(ipx) - 1
649	ila_pcx(ipx+1) = ila_pcx(ipx+1) + 1
650	ELSE
651	glo_rank(j+1) = glo_rank(j)
652	ila_pcx(ipx) = ila_pcx(ipx) + 1
653	ila_pcx(ipx+1) = ila_pcx(ipx+1) - 1
654	END IF
655	j = j + 1
656	END IF
657	j = j + 1
658	IF (ila_pcx(ipx) == 0) EXIT
659	END DO
660	END DO
661	END DO
662	! Move to next npy slice
663	ib = ie
664	END DO
665
666	! Count local elements and assign local values
667	nbp = COUNT(glo_rank == mpi_rank)
668	return_ni=nbp
669	return_nj=1
670
671	ALLOCATE(return_lon(0:nbp-1))
672	ALLOCATE(return_lat(0:nbp-1))
673	ALLOCATE(return_mask(0:nbp-1))
674	ALLOCATE(return_index(0:nbp-1))
675	ALLOCATE(i_index(0:nbp-1))
676	ALLOCATE(bounds_lon(ncrn,0:nbp-1))
677	ALLOCATE(bounds_lat(ncrn,0:nbp-1))
678	ALLOCATE(srf_loc2d(0:0,0:nbp-1))
679
680	j = 0
681	DO i = 0, ni_glo-1
682	IF (glo_rank(i) == mpi_rank) THEN
683	return_lon(j) = lon_glo2D(i,0)
684	return_lat(j) = lat_glo2D(i,0)
685	return_mask(j) = msk_glo2D(i,0) == 0
686	return_index(j) = j
687	i_index(j) = i
688	DO ic = 1, ncrn
689	bounds_lon(ic,j) = clo_glo2d(i,0,ic)
690	bounds_lat(ic,j) = cla_glo2d(i,0,ic)
691	END DO
692	!AP YANN Notice the inversion of indexes for 1d areas stored in 2d
693	srf_loc2d(0,j) = srf_glo2d(i,0)
694	j = j + 1 ! Next local cell
695	END IF
696	END DO
697
698	IF ( .NOT. params%domain_mask) return_mask(:) = .TRUE.
699
700	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
701	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="unstructured", data_dim=1)
702	CALL xios_set_domain_attr(TRIM(domain_id), ni_glo=ni_glo, ibegin=0, ni=nbp, i_index=i_index)
703	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_1d=return_lon, latvalue_1d=return_lat, mask_1d=return_mask)
704	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat, nvertex=ncrn)
705	CALL xios_set_domain_attr(TRIM(domain_id), area=srf_loc2d)
706	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
707	ENDIF
708
709	DEALLOCATE(glo_rank, ila_pcx, i_index, srf_loc2d, bounds_lon, bounds_lat)
710
711	END SUBROUTINE init_domain_ssea
712
713	SUBROUTINE init_domain_lmdz(domain_id, comm, params, return_ni, return_nj, &
714	return_lon,return_lat,return_mask, return_index)
715	IMPLICIT NONE
716	CHARACTER(LEN=*) :: domain_id
717	TYPE(tmodel_params) :: params
718	INTEGER :: comm
719	INTEGER :: return_ni
720	INTEGER :: return_nj
721	DOUBLE PRECISION, POINTER :: return_lon(:)
722	DOUBLE PRECISION, POINTER :: return_lat(:)
723	LOGICAL, POINTER :: return_mask(:)
724	INTEGER, POINTER :: return_index(:)
725	INTEGER :: domain_proc_rank, domain_proc_size
726	INTEGER :: axis_proc_rank, axis_proc_size
727
728	INTEGER :: mpi_rank, mpi_size
729	INTEGER :: ierr
730	INTEGER :: ni,nj,ni_glo,nj_glo
731	INTEGER :: ibegin,jbegin
732	INTEGER :: nbp,nbp_glo, offset
733	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), lon(:), lat(:)
734	LOGICAL,ALLOCATABLE :: mask(:)
735	LOGICAL,ALLOCATABLE :: dom_mask(:)
736	INTEGER :: i,j
737
738	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
739	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
740
741	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank)
742
743	ni_glo=params%ni
744	nj_glo=params%nj
745	nbp_glo=ni_glo*nj_glo
746	nbp=nbp_glo/domain_proc_size
747	IF (domain_proc_rank < MOD(nbp_glo,domain_proc_size) ) THEN
748	nbp=nbp+1
749	offset=nbp*domain_proc_rank
750	ELSE
751	offset=nbp*domain_proc_rank + MOD(nbp_glo, domain_proc_size)
752	ENDIF
753
754
755	ibegin=0 ; ni=ni_glo
756	jbegin=offset / ni_glo
757
758	nj = (offset+nbp-1) / ni_glo - jbegin + 1
759
760	offset=offset-jbegin*ni
761
762	ALLOCATE(lon(0:ni-1), lat(0:nj-1), mask(0:ninj-1), dom_mask(0:ninj-1))
763	mask(:)=.FALSE.
764	mask(offset:offset+nbp-1)=.TRUE.
765
766	ALLOCATE(lon_glo(0:ni_glo-1), lat_glo(0:nj_glo-1))
767
768	DO i=0,ni_glo-1
769	lon_glo(i)=-180+(i+0.5)*(360./ni_glo)
770	ENDDO
771
772	DO j=0,nj_glo-1
773	lat_glo(j)=-90+(j+0.5)*(180./nj_glo)
774	ENDDO
775
776	lon(:)=lon_glo(:)
777	lat(:)=lat_glo(jbegin:jbegin+nj-1)
778
779	ALLOCATE(return_lon(0:ni*nj-1))
780	ALLOCATE(return_lat(0:ni*nj-1))
781	ALLOCATE(return_mask(0:ni*nj-1))
782	ALLOCATE(return_index(0:ni*nj-1))
783
784	DO j=0,nj-1
785	DO i=0,ni-1
786	return_lon(i+j*ni)=lon(i)
787	return_lat(i+j*ni)=lat(j)
788	return_index(i+jni)=i+jni
789	ENDDO
790	ENDDO
791
792	CALL set_domain_mask(params, return_lon,return_lat, dom_mask)
793
794	return_mask = mask .AND. dom_mask
795
796	return_ni=ni
797	return_nj=nj
798
799	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
800	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="rectilinear", ni_glo=ni_glo, ibegin=ibegin, ni=ni, nj_glo=nj_glo, &
801	jbegin=jbegin, nj=nj)
802	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=2, lonvalue_1d=lon, latvalue_1d=lat, mask_1d=return_mask)
803	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
804	ENDIF
805
806	END SUBROUTINE init_domain_lmdz
807
808	SUBROUTINE init_domain_bggd(domain_id, comm, params, return_ni, return_nj, &
809	return_lon,return_lat,return_mask, return_index)
810	IMPLICIT NONE
811	CHARACTER(LEN=*) :: domain_id
812	TYPE(tmodel_params) :: params
813	INTEGER :: comm
814	INTEGER :: return_ni
815	INTEGER :: return_nj
816	DOUBLE PRECISION, POINTER :: return_lon(:)
817	DOUBLE PRECISION, POINTER :: return_lat(:)
818	LOGICAL, POINTER :: return_mask(:)
819	INTEGER, POINTER :: return_index(:)
820
821	INTEGER :: mpi_rank, mpi_size
822	INTEGER :: ierr
823
824	DOUBLE PRECISION, POINTER :: lon_glo2d(:,:), lat_glo2d(:,:)
825	DOUBLE PRECISION, POINTER :: clo_glo2d(:,:,:), cla_glo2d(:,:,:)
826	INTEGER, POINTER :: msk_glo2d(:,:)
827	DOUBLE PRECISION, POINTER :: srf_glo2d(:,:)
828
829	INTEGER :: ni,nj,ni_glo,nj_glo, ncrn
830	INTEGER :: ibegin,jbegin
831	INTEGER :: nbp,nbp_glo, offset
832	DOUBLE PRECISION, ALLOCATABLE :: lon2d(:,:), lat2d(:,:), srf_loc2d(:,:)
833	DOUBLE PRECISION, ALLOCATABLE :: bounds_lon(:,:,:), bounds_lat(:,:,:)
834	LOGICAL,ALLOCATABLE :: mask(:)
835	LOGICAL,ALLOCATABLE :: dom_mask(:)
836	INTEGER :: i,j,ij,ic
837
838	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
839	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
840
841	CALL read_oasis_grid(params, ni_glo, nj_glo, ncrn,&
842	& lon_glo2d, lat_glo2d, clo_glo2d, cla_glo2d, msk_glo2d, srf_glo2d)
843
844	nbp_glo=ni_glo*nj_glo
845	nbp=nbp_glo/mpi_size
846	IF (mpi_rank < MOD(nbp_glo,mpi_size) ) THEN
847	nbp=nbp+1
848	offset=nbp*mpi_rank
849	ELSE
850	offset=nbp*mpi_rank + MOD(nbp_glo, mpi_size)
851	ENDIF
852
853	ibegin=0 ; ni=ni_glo
854	jbegin=offset / ni_glo
855
856	nj = (offset+nbp-1) / ni_glo - jbegin + 1
857
858	offset=offset-jbegin*ni
859
860	ALLOCATE(mask(0:ninj-1), dom_mask(0:ninj-1))
861	ALLOCATE(lon2d(0:ni-1,0:nj-1), lat2d(0:ni-1,0:nj-1))
862	mask(:)=.FALSE.
863	mask(offset:offset+nbp-1)=.TRUE.
864
865	lon2d(:,:)=lon_glo2d(ibegin:ibegin+ni-1,jbegin:jbegin+nj-1)
866	lat2d(:,:)=lat_glo2d(ibegin:ibegin+ni-1,jbegin:jbegin+nj-1)
867
868	ALLOCATE(bounds_lon(ncrn,0:ni-1,0:nj-1))
869	ALLOCATE(bounds_lat(ncrn,0:ni-1,0:nj-1))
870	ALLOCATE(srf_loc2D(0:ni-1,0:nj-1))
871
872	srf_loc2D(:,:)=srf_glo2d(ibegin:ibegin+ni-1,jbegin:jbegin+nj-1)
873
874	ALLOCATE(return_lon(0:ni*nj-1))
875	ALLOCATE(return_lat(0:ni*nj-1))
876	ALLOCATE(return_mask(0:ni*nj-1))
877	ALLOCATE(return_index(0:ni*nj-1))
878
879	DO j=0,nj-1
880	DO i=0,ni-1
881	ij = i+j*ni
882	return_lon(ij)=lon_glo2d(ibegin+i,jbegin+j)
883	return_lat(ij)=lat_glo2d(ibegin+i,jbegin+j)
884	dom_mask(ij) = msk_glo2d(ibegin+i,jbegin+j) == 0
885	return_index(ij)=ij
886	DO ic = 1, ncrn
887	bounds_lon(ic,i,j) = clo_glo2d(ibegin+i,jbegin+j,ic)
888	bounds_lat(ic,i,j) = cla_glo2d(ibegin+i,jbegin+j,ic)
889	END DO
890	ENDDO
891	ENDDO
892
893	return_mask = mask .AND. dom_mask
894	IF ( .NOT. params%domain_mask ) return_mask = mask
895
896	return_ni=ni
897	return_nj=nj
898
899	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
900	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="rectilinear", data_dim=2)
901	CALL xios_set_domain_attr(TRIM(domain_id), ni_glo=ni_glo, ibegin=ibegin, ni=ni)
902	CALL xios_set_domain_attr(TRIM(domain_id), nj_glo=nj_glo, jbegin=jbegin, nj=nj)
903	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_2d=lon2d, latvalue_2d=lat2d, mask_1d=return_mask)
904	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat, nvertex=ncrn)
905	CALL xios_set_domain_attr(TRIM(domain_id), area=srf_loc2d)
906	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
907	ENDIF
908
909	DEALLOCATE(lon2d, lat2d, srf_loc2d, bounds_lon, bounds_lat, mask, dom_mask)
910
911	END SUBROUTINE init_domain_bggd
912
913
914	SUBROUTINE init_domain_nemo(domain_id, comm, params, return_ni, return_nj, &
915	return_lon,return_lat,return_mask, return_index)
916	IMPLICIT NONE
917	CHARACTER(LEN=*) :: domain_id
918	TYPE(tmodel_params) :: params
919	INTEGER :: comm
920	INTEGER :: return_ni
921	INTEGER :: return_nj
922	DOUBLE PRECISION, POINTER :: return_lon(:)
923	DOUBLE PRECISION, POINTER :: return_lat(:)
924	LOGICAL, POINTER :: return_mask(:)
925	INTEGER, POINTER :: return_index(:)
926	INTEGER :: domain_proc_rank, domain_proc_size
927	INTEGER :: axis_proc_rank, axis_proc_size
928
929	INTEGER :: mpi_rank, mpi_size
930	INTEGER :: ierr
931	INTEGER :: ni,nj,ni_glo,nj_glo
932	INTEGER :: ibegin,jbegin
933	INTEGER :: offset_i, offset_j
934	DOUBLE PRECISION, ALLOCATABLE :: lon_glo(:), lat_glo(:), bounds_lon_glo(:,:), bounds_lat_glo(:,:)
935	DOUBLE PRECISION, ALLOCATABLE :: lon(:,:), lat(:,:), bounds_lon(:,:,:), bounds_lat(:,:,:)
936	INTEGER :: i,j, ij, n, rank
937	INTEGER :: nproc_i, nproc_j, nholes, pos_hole
938	INTEGER,ALLOCATABLE :: size_i(:), begin_i(:), size_j(:), begin_j(:)
939
940	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
941	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
942
943	CALL get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank)
944	ni_glo=params%ni
945	nj_glo=params%nj
946
947	n=INT(SQRT(mpi_size*1.))
948	nproc_i=n
949	nproc_j=n
950	IF ( n*n == mpi_size) THEN
951	! do nothing
952	ELSE IF ( (n+1)*n < mpi_size) THEN
953	nproc_i=nproc_i+1
954	nproc_j=nproc_j+1
955	ELSE
956	nproc_i=nproc_i+1
957	ENDIF
958
959	nholes=nproc_i*nproc_j-mpi_size
960
961	ALLOCATE(size_i(0:nproc_i-1))
962	ALLOCATE(begin_i(0:nproc_i-1))
963	DO i=0,nproc_i-1
964	size_i(i)=ni_glo/nproc_i
965	IF (i<MOD(ni_glo,nproc_i)) size_i(i)=size_i(i)+1
966	IF (i==0) THEN
967	begin_i(i)=0
968	ELSE
969	begin_i(i)=begin_i(i-1)+size_i(i-1)
970	ENDIF
971	ENDDO
972
973	ALLOCATE(size_j(0:nproc_j-1))
974	ALLOCATE(begin_j(0:nproc_j-1))
975	DO j=0,nproc_j-1
976	size_j(j)=nj_glo/nproc_j
977	IF (j<MOD(nj_glo,nproc_j)) size_j(j)=size_j(j)+1
978	IF (j==0) THEN
979	begin_j(j)=0
980	ELSE
981	begin_j(j)=begin_j(j-1)+size_j(j-1)
982	ENDIF
983	ENDDO
984
985
986	pos_hole=0
987	rank=0
988	DO j=0,nproc_j-1
989	DO i=0,nproc_i-1
990
991	ij = j*nproc_i + i
992	IF (pos_hole<nholes) THEN
993	IF ( MOD(ij,(nproc_i*nproc_j/nholes)) == 0) THEN
994	pos_hole=pos_hole+1
995	CYCLE
996	ENDIF
997	ENDIF
998
999	IF (mpi_rank==rank) THEN
1000	ibegin = begin_i(i)
1001	ni = size_i(i)
1002	jbegin = begin_j(j)
1003	nj = size_j(j)
1004	ENDIF
1005	rank=rank+1
1006	ENDDO
1007	ENDDO
1008
1009	ALLOCATE(lon_glo(0:ni_glo-1), lat_glo(0:nj_glo-1))
1010	ALLOCATE(bounds_lon_glo(4,0:ni_glo-1), bounds_lat_glo(4,0:nj_glo-1))
1011
1012	DO i=0,ni_glo-1
1013	lon_glo(i)=-180+(i+0.5)*(360./ni_glo)
1014	bounds_lon_glo(1,i)= -180+(i)*(360./ni_glo)
1015	bounds_lon_glo(2,i)= -180+(i+1)*(360./ni_glo)
1016	ENDDO
1017
1018	DO j=0,nj_glo-1
1019	lat_glo(j)=-90+(j+0.5)*(180./nj_glo)
1020	bounds_lat_glo(1,j)= -90+(j)*(180./nj_glo)
1021	bounds_lat_glo(2,j)= -90+(j+1)*(180./nj_glo)
1022	ENDDO
1023
1024	offset_i=2 ! halo of 2 on i
1025	offset_j=1 ! halo of 1 on j
1026
1027	ALLOCATE(lon(0:ni-1,0:nj-1))
1028	ALLOCATE(lat(0:ni-1,0:nj-1))
1029	ALLOCATE(bounds_lon(4,0:ni-1,0:nj-1))
1030	ALLOCATE(bounds_lat(4,0:ni-1,0:nj-1))
1031
1032	ALLOCATE(return_lon(0:ni*nj-1))
1033	ALLOCATE(return_lat(0:ni*nj-1))
1034	ALLOCATE(return_mask(0:ni*nj-1))
1035	ALLOCATE(return_index(0:(ni+2offset_i)(nj+2*offset_j)-1))
1036
1037	return_index=-1
1038	DO j=0,nj-1
1039	DO i=0,ni-1
1040	ij=j*ni+i
1041	return_lon(ij)=lon_glo(ibegin+i)
1042	return_lat(ij)=lat_glo(jbegin+j)
1043	lon(i,j)=return_lon(ij)
1044	lat(i,j)=return_lat(ij)
1045	bounds_lon(1,i,j)=bounds_lon_glo(2,ibegin+i)
1046	bounds_lon(2,i,j)=bounds_lon_glo(1,ibegin+i)
1047	bounds_lon(3,i,j)=bounds_lon_glo(1,ibegin+i)
1048	bounds_lon(4,i,j)=bounds_lon_glo(2,ibegin+i)
1049	bounds_lat(1,i,j)=bounds_lat_glo(1,jbegin+j)
1050	bounds_lat(2,i,j)=bounds_lat_glo(1,jbegin+j)
1051	bounds_lat(3,i,j)=bounds_lat_glo(2,jbegin+j)
1052	bounds_lat(4,i,j)=bounds_lat_glo(2,jbegin+j)
1053
1054	ij=(j+offset_j)(ni+2offset_i)+i+offset_i
1055	return_index(ij)=i+j*ni
1056	ENDDO
1057	ENDDO
1058
1059	CALL set_domain_mask(params, return_lon,return_lat, return_mask)
1060
1061	return_ni=ni
1062	return_nj=nj
1063
1064	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1065	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="curvilinear", data_dim=2)
1066	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)
1067	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)
1068	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=2, lonvalue_2d=lon, latvalue_2d=lat, mask_1d=return_mask)
1069	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat, nvertex=4)
1070	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
1071	ENDIF
1072
1073	END SUBROUTINE init_domain_nemo
1074
1075
1076	SUBROUTINE init_domain_orca(domain_id, comm, params, return_ni, return_nj, &
1077	return_lon,return_lat,return_mask, return_index)
1078	IMPLICIT NONE
1079	CHARACTER(LEN=*) :: domain_id
1080	TYPE(tmodel_params) :: params
1081	INTEGER :: comm
1082	INTEGER :: return_ni
1083	INTEGER :: return_nj
1084	DOUBLE PRECISION, POINTER :: return_lon(:)
1085	DOUBLE PRECISION, POINTER :: return_lat(:)
1086	LOGICAL, POINTER :: return_mask(:)
1087	INTEGER, POINTER :: return_index(:)
1088	INTEGER :: mpi_rank, mpi_size
1089	INTEGER :: ierr
1090
1091	DOUBLE PRECISION, POINTER :: lon_glo2d(:,:), lat_glo2d(:,:)
1092	DOUBLE PRECISION, POINTER :: clo_glo2d(:,:,:), cla_glo2d(:,:,:)
1093	INTEGER, POINTER :: msk_glo2d(:,:)
1094	DOUBLE PRECISION, POINTER :: srf_glo2d(:,:)
1095
1096	INTEGER :: ni,nj,ni_glo,nj_glo, ncrn
1097	INTEGER :: ibegin,jbegin
1098	INTEGER :: offset_i, offset_j
1099	DOUBLE PRECISION, ALLOCATABLE :: lon(:,:), lat(:,:), bounds_lon(:,:,:), bounds_lat(:,:,:), srf_loc2d(:,:)
1100	INTEGER :: i, j, ij, rank
1101	INTEGER :: nproc_i, nproc_j
1102	INTEGER,ALLOCATABLE :: size_i(:), begin_i(:), size_j(:), begin_j(:)
1103
1104	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1105	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1106
1107	CALL read_oasis_grid(params, ni_glo, nj_glo, ncrn,&
1108	& lon_glo2d, lat_glo2d, clo_glo2d, cla_glo2d, msk_glo2d, srf_glo2d)
1109
1110	nproc_j=INT(SQRT(mpi_size*1.))
1111	DO WHILE(MOD(mpi_size,nproc_j) /= 0)
1112	nproc_j = nproc_j - 1
1113	END DO
1114	nproc_i = mpi_size / nproc_j
1115
1116	ALLOCATE(size_i(0:nproc_i-1))
1117	ALLOCATE(begin_i(0:nproc_i-1))
1118	DO i=0,nproc_i-1
1119	size_i(i)=ni_glo/nproc_i
1120	IF (i<MOD(ni_glo,nproc_i)) size_i(i)=size_i(i)+1
1121	IF (i==0) THEN
1122	begin_i(i)=0
1123	ELSE
1124	begin_i(i)=begin_i(i-1)+size_i(i-1)
1125	ENDIF
1126	ENDDO
1127
1128	ALLOCATE(size_j(0:nproc_j-1))
1129	ALLOCATE(begin_j(0:nproc_j-1))
1130	DO j=0,nproc_j-1
1131	size_j(j)=nj_glo/nproc_j
1132	IF (j<MOD(nj_glo,nproc_j)) size_j(j)=size_j(j)+1
1133	IF (j==0) THEN
1134	begin_j(j)=0
1135	ELSE
1136	begin_j(j)=begin_j(j-1)+size_j(j-1)
1137	ENDIF
1138	ENDDO
1139
1140	rank=0
1141	DO j=0,nproc_j-1
1142	DO i=0,nproc_i-1
1143	IF (mpi_rank==rank) THEN
1144	ibegin = begin_i(i)
1145	ni = size_i(i)
1146	jbegin = begin_j(j)
1147	nj = size_j(j)
1148	ENDIF
1149	rank=rank+1
1150	ENDDO
1151	ENDDO
1152
1153	offset_i=2 ! halo of 2 on i
1154	offset_j=1 ! halo of 1 on j
1155
1156	ALLOCATE(lon(0:ni-1,0:nj-1))
1157	ALLOCATE(lat(0:ni-1,0:nj-1))
1158	ALLOCATE(bounds_lon(ncrn,0:ni-1,0:nj-1))
1159	ALLOCATE(bounds_lat(ncrn,0:ni-1,0:nj-1))
1160	ALLOCATE(srf_loc2D(0:ni-1,0:nj-1))
1161
1162	ALLOCATE(return_lon(0:ni*nj-1))
1163	ALLOCATE(return_lat(0:ni*nj-1))
1164	ALLOCATE(return_mask(0:ni*nj-1))
1165	ALLOCATE(return_index(0:(ni+2offset_i)(nj+2*offset_j)-1))
1166
1167
1168	return_index=-1
1169	DO j=0,nj-1
1170	DO i=0,ni-1
1171	ij=j*ni+i
1172	return_lon(ij)=lon_glo2d(ibegin+i,jbegin+j)
1173	return_lat(ij)=lat_glo2d(ibegin+i,jbegin+j)
1174	return_mask(ij) = msk_glo2d(ibegin+i,jbegin+j) == 0
1175	srf_loc2d(i,j) = srf_glo2d(ibegin+i,jbegin+j)
1176	lon(i,j)=return_lon(ij)
1177	lat(i,j)=return_lat(ij)
1178	DO ij = 1, ncrn
1179	bounds_lon(ij,i,j) = clo_glo2d(ibegin+i,jbegin+j,ij)
1180	bounds_lat(ij,i,j) = cla_glo2d(ibegin+i,jbegin+j,ij)
1181	END DO
1182
1183	ij=(j+offset_j)(ni+2offset_i)+i+offset_i
1184	return_index(ij)=i+j*ni
1185	ENDDO
1186	ENDDO
1187
1188	return_ni=ni
1189	return_nj=nj
1190
1191	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1192	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="curvilinear", data_dim=2)
1193	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)
1194	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)
1195	CALL xios_set_domain_attr(TRIM(domain_id), data_dim=2, lonvalue_2d=lon, latvalue_2d=lat, mask_1d=return_mask)
1196	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat, nvertex=ncrn)
1197	CALL xios_set_domain_attr(TRIM(domain_id), area=srf_loc2d)
1198	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
1199	ENDIF
1200
1201	DEALLOCATE(lon, lat, bounds_lon, bounds_lat, srf_loc2d, size_i, begin_i, size_j, begin_j)
1202
1203	END SUBROUTINE init_domain_orca
1204
1205	SUBROUTINE init_domain_icos(domain_id, comm, params, return_ni, return_nj, &
1206	return_lon,return_lat,return_mask, return_index)
1207	IMPLICIT NONE
1208	CHARACTER(LEN=*) :: domain_id
1209	TYPE(tmodel_params) :: params
1210	INTEGER :: comm
1211	INTEGER :: return_ni
1212	INTEGER :: return_nj
1213	DOUBLE PRECISION, POINTER :: return_lon(:)
1214	DOUBLE PRECISION, POINTER :: return_lat(:)
1215	LOGICAL, POINTER :: return_mask(:)
1216	INTEGER, POINTER :: return_index(:)
1217	INTEGER :: mpi_rank, mpi_size
1218	INTEGER :: ierr
1219
1220	DOUBLE PRECISION, POINTER :: lon_glo2d(:,:), lat_glo2d(:,:)
1221	DOUBLE PRECISION, POINTER :: clo_glo2d(:,:,:), cla_glo2d(:,:,:)
1222	INTEGER, POINTER :: msk_glo2d(:,:)
1223	DOUBLE PRECISION, POINTER :: srf_glo2d(:,:)
1224
1225	INTEGER :: ni_glo,nj_glo, ncrn
1226	INTEGER :: nbp,nbp_glo, offset
1227
1228	DOUBLE PRECISION, ALLOCATABLE :: lon(:), lat(:), srf_loc2d(:,:)
1229	DOUBLE PRECISION, ALLOCATABLE :: bounds_lon(:,:), bounds_lat(:,:)
1230
1231	INTEGER :: i,ic
1232
1233	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1234	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1235
1236	CALL read_oasis_grid(params, ni_glo, nj_glo, ncrn,&
1237	& lon_glo2d, lat_glo2d, clo_glo2d, cla_glo2d, msk_glo2d, srf_glo2d)
1238
1239
1240	nbp_glo=ni_glo
1241	nbp=nbp_glo/mpi_size
1242	IF (mpi_rank < MOD(nbp_glo,mpi_size) ) THEN
1243	nbp=nbp+1
1244	offset=nbp*mpi_rank
1245	ELSE
1246	offset=nbp*mpi_rank + MOD(nbp_glo, mpi_size)
1247	ENDIF
1248
1249	ALLOCATE(lat(0:nbp-1))
1250	ALLOCATE(lon(0:nbp-1))
1251	ALLOCATE(bounds_lon(ncrn,0:nbp-1))
1252	ALLOCATE(bounds_lat(ncrn,0:nbp-1))
1253	ALLOCATE(return_lon(0:nbp-1))
1254	ALLOCATE(return_lat(0:nbp-1))
1255	ALLOCATE(return_index(0:nbp-1))
1256	ALLOCATE(return_mask(0:nbp-1))
1257	ALLOCATE(srf_loc2d(0:0,0:nbp-1))
1258
1259	DO i=0,nbp-1
1260	lat(i)=lat_glo2d(i+offset,0)
1261	lon(i)=lon_glo2d(i+offset,0)
1262	DO ic = 1, ncrn
1263	bounds_lon(ic,i) = clo_glo2d(i+offset,0,ic)
1264	bounds_lat(ic,i) = cla_glo2d(i+offset,0,ic)
1265	END DO
1266	return_mask(i) = msk_glo2d(i+offset,0) == 0
1267	!AP YANN Notice the inversion of indexes for 1d areas stored in 2d
1268	srf_loc2d(0,i) = srf_glo2d(i+offset,0)
1269	return_index(i)=i
1270	ENDDO
1271
1272	IF ( .NOT. params%domain_mask ) return_mask(:) = .TRUE.
1273
1274	return_lon=lon
1275	return_lat=lat
1276
1277	return_ni=nbp
1278	return_nj=1
1279
1280	IF (xios_is_valid_domain(TRIM(domain_id))) THEN
1281	CALL xios_set_domain_attr(TRIM(domain_id), TYPE="unstructured", data_dim=1)
1282	CALL xios_set_domain_attr(TRIM(domain_id), ni_glo=ni_glo, ibegin=offset, ni=nbp)
1283	CALL xios_set_domain_attr(TRIM(domain_id), lonvalue_1d=lon, latvalue_1d=lat, mask_1d=return_mask)
1284	CALL xios_set_domain_attr(TRIM(domain_id), bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat, nvertex=ncrn)
1285	CALL xios_set_domain_attr(TRIM(domain_id), area=srf_loc2d)
1286	CALL xios_set_domain_attr(TRIM(domain_id), radius=1.0_8)
1287	ENDIF
1288
1289	DEALLOCATE(lon, lat, srf_loc2d, bounds_lon, bounds_lat)
1290
1291	END SUBROUTINE init_domain_icos
1292
1293
1294	SUBROUTINE set_domain_mask(params, lon,lat, mask)
1295	IMPLICIT NONE
1296	TYPE(tmodel_params) :: params
1297	DOUBLE PRECISION :: lon(:)
1298	DOUBLE PRECISION :: lat(:)
1299	LOGICAL :: mask(:)
1300
1301	mask(:)=.TRUE.
1302	IF (params%domain_mask) THEN
1303	WHERE (lon(:)-2lat(:)>-10 .AND. lon(:)-2lat(:) <10) mask(:)=.FALSE.
1304	WHERE (2lat(:)+lon(:)>-10 .AND. 2lat(:)+lon(:)<10) mask(:)=.FALSE.
1305	ENDIF
1306
1307	END SUBROUTINE set_domain_mask
1308
1309
1310	SUBROUTINE init_field2D_vortex(params, lon, lat, mask, return_field)
1311	IMPLICIT NONE
1312	TYPE(tmodel_params) :: params
1313	DOUBLE PRECISION, POINTER :: lon(:)
1314	DOUBLE PRECISION, POINTER :: lat(:)
1315	LOGICAL, POINTER :: mask(:)
1316	DOUBLE PRECISION, POINTER :: return_field(:)
1317
1318	DOUBLE PRECISION, PARAMETER :: dp_pi=3.14159265359
1319	DOUBLE PRECISION, PARAMETER :: dLon0 = 5.5
1320	DOUBLE PRECISION, PARAMETER :: dLat0 = 0.2
1321	DOUBLE PRECISION, PARAMETER :: dR0 = 3.0
1322	DOUBLE PRECISION, PARAMETER :: dD = 5.0
1323	DOUBLE PRECISION, PARAMETER :: dT = 6.0
1324	DOUBLE PRECISION :: dp_length, dp_conv
1325	INTEGER :: i,xy
1326
1327	DOUBLE PRECISION :: dSinC, dCosC, dCosT, dSinT
1328	DOUBLE PRECISION :: dTrm, dX, dY, dZ
1329	DOUBLE PRECISION :: dlon, dlat
1330	DOUBLE PRECISION :: dRho, dVt, dOmega
1331
1332	xy=SIZE(mask)
1333
1334	ALLOCATE(return_field(0:xy-1))
1335
1336	dp_conv = dp_pi/180.
1337	dSinC = SIN( dLat0 )
1338	dCosC = COS( dLat0 )
1339
1340	DO i=0,xy-1
1341	IF (mask(i)) THEN
1342
1343	! Find the rotated longitude and latitude of a point on a sphere
1344	! with pole at (dLon0, dLat0).
1345	dCosT = COS( lat(i)*dp_conv )
1346	dSinT = SIN( lat(i)*dp_conv )
1347
1348	dTrm = dCosT * COS( lon(i)*dp_conv - dLon0 )
1349	dX = dSinC * dTrm - dCosC * dSinT
1350	dY = dCosT * SIN( lon(i)*dp_conv - dLon0 )
1351	dZ = dSinC * dSinT + dCosC * dTrm
1352
1353	dlon = ATAN2( dY, dX )
1354	IF( dlon < 0.0 ) dlon = dlon + 2.0 * dp_pi
1355	dlat = ASIN( dZ )
1356
1357	dRho = dR0 * COS(dlat)
1358	dVt = 3.0 * SQRT(3.0)/2.0/COSH(dRho)/COSH(dRho)*TANH(dRHo)
1359	IF (dRho == 0.0) THEN
1360	dOmega = 0.0
1361	ELSE
1362	dOmega = dVt / dRho
1363	END IF
1364
1365	RETURN_FIELD(i) = 2.0 * ( 1.0 + TANH( dRho / dD * SIN( dLon - dOmega * dT ) ) )
1366
1367	END IF
1368	END DO
1369
1370	END SUBROUTINE init_field2D_vortex
1371
1372	SUBROUTINE init_field2D_academic(params, lon, lat, mask, return_field)
1373	IMPLICIT NONE
1374	TYPE(tmodel_params) :: params
1375	DOUBLE PRECISION, POINTER :: lon(:)
1376	DOUBLE PRECISION, POINTER :: lat(:)
1377	LOGICAL, POINTER :: mask(:)
1378	DOUBLE PRECISION, POINTER :: return_field(:)
1379
1380	DOUBLE PRECISION, PARAMETER :: coef=2., dp_pi=3.14159265359
1381	DOUBLE PRECISION :: dp_length, dp_conv
1382	INTEGER :: i,xy
1383
1384	! Analytical Gulf Stream
1385	DOUBLE PRECISION :: gf_coef, gf_ori_lon, gf_ori_lat, &
1386	& gf_end_lon, gf_end_lat, gf_dmp_lon, gf_dmp_lat
1387	DOUBLE PRECISION :: gf_per_lon
1388	DOUBLE PRECISION :: dx, dy, dr, dth, dc, dr0, dr1
1389
1390	! Parameters for analytical function
1391	dp_length= 1.2*dp_pi
1392	dp_conv=dp_pi/180.
1393	gf_coef = params%gf_coef ! Coefficient for Gulf Stream term (0.0 = no Gulf Stream)
1394	gf_ori_lon = -80.0 ! Origin of the Gulf Stream (longitude in deg)
1395	gf_ori_lat = 25.0 ! Origin of the Gulf Stream (latitude in deg)
1396	gf_end_lon = -1.8 ! End point of the Gulf Stream (longitude in deg)
1397	gf_end_lat = 50.0 ! End point of the Gulf Stream (latitude in deg)
1398	gf_dmp_lon = -25.5 ! Point of the Gulf Stream decrease (longitude in deg)
1399	gf_dmp_lat = 55.5 ! Point of the Gulf Stream decrease (latitude in deg)
1400
1401	xy=SIZE(mask)
1402
1403	ALLOCATE(return_field(0:xy-1))
1404
1405	dr0 = SQRT(((gf_end_lon - gf_ori_lon)dp_conv)*2 + &
1406	& ((gf_end_lat - gf_ori_lat)dp_conv)*2)
1407	dr1 = SQRT(((gf_dmp_lon - gf_ori_lon)dp_conv)*2 + &
1408	& ((gf_dmp_lat - gf_ori_lat)dp_conv)*2)
1409
1410	DO i=0,xy-1
1411	IF (mask(i)) THEN
1412	return_field(i)=(coef-COS(dp_pi(ACOS(COS(lat(i)dp_conv)*&
1413	COS(lon(i)*dp_conv))/dp_length)))
1414	gf_per_lon = lon(i)
1415	IF (gf_per_lon > 180.0) gf_per_lon = gf_per_lon - 360.0
1416	IF (gf_per_lon < -180.0) gf_per_lon = gf_per_lon + 360.0
1417	dx = (gf_per_lon - gf_ori_lon)*dp_conv
1418	dy = (lat(i) - gf_ori_lat)*dp_conv
1419	dr = SQRT(dxdx + dydy)
1420	dth = ATAN2(dy, dx)
1421	dc = 1.3 * gf_coef
1422	IF (dr > dr0) dc = 0.0
1423	IF (dr > dr1) dc = dc * COS(dp_pi0.5(dr-dr1)/(dr0-dr1))
1424	return_field(i) = return_field(i) + &
1425	& (MAX(1000.0SIN(0.4(0.5*dr+dth) + &
1426	& 0.007COS(50.0dth) + 0.37dp_pi),999.0) - 999.0) &
1427	& dc
1428	ENDIF
1429	ENDDO
1430
1431	END SUBROUTINE init_field2D_academic
1432
1433	SUBROUTINE init_field2D_one(mask, return_field)
1434
1435	IMPLICIT NONE
1436	LOGICAL, POINTER :: mask(:)
1437	DOUBLE PRECISION, POINTER :: return_field(:)
1438
1439	INTEGER :: xy
1440
1441	xy=SIZE(mask)
1442
1443	ALLOCATE(return_field(0:xy-1))
1444	return_field(:)=1.0
1445
1446	END SUBROUTINE init_field2D_one
1447
1448	SUBROUTINE init_field2D_rank(comm, mask, return_field)
1449
1450	IMPLICIT NONE
1451	INTEGER :: comm
1452	LOGICAL, POINTER :: mask(:)
1453	DOUBLE PRECISION, POINTER :: return_field(:)
1454
1455	INTEGER :: rank,ierr
1456	INTEGER :: xy
1457
1458	CALL MPI_COMM_RANK(comm,rank,ierr)
1459
1460	xy=SIZE(mask)
1461
1462	ALLOCATE(return_field(0:xy-1))
1463	return_field(:)=rank
1464
1465	END SUBROUTINE init_field2D_rank
1466
1467	SUBROUTINE init_field2D_mask(comm, mask, return_field)
1468
1469	IMPLICIT NONE
1470	INTEGER :: comm
1471	LOGICAL, POINTER :: mask(:)
1472	DOUBLE PRECISION, POINTER :: return_field(:)
1473
1474	INTEGER :: rank,ierr
1475	INTEGER :: xy
1476
1477	CALL MPI_COMM_RANK(comm,rank,ierr)
1478
1479	xy=SIZE(mask)
1480
1481	ALLOCATE(return_field(0:xy-1))
1482	WHERE (mask)
1483	return_field(:) = 1.0
1484	ELSEWHERE
1485	return_field(:) = 0.0
1486	END WHERE
1487
1488	END SUBROUTINE init_field2D_mask
1489
1490	SUBROUTINE read_oasis_grid(params, ni_glo, nj_glo, ncrn,&
1491	& lon_glo2d, lat_glo2d, clo_glo2d, cla_glo2d, msk_glo2d, srf_glo2d)
1492	IMPLICIT NONE
1493	TYPE(tmodel_params), INTENT(IN) :: params
1494	INTEGER, INTENT(OUT) :: ni_glo, nj_glo, ncrn
1495	DOUBLE PRECISION, POINTER :: lon_glo2d(:,:), lat_glo2d(:,:)
1496	DOUBLE PRECISION, POINTER :: clo_glo2d(:,:,:), cla_glo2d(:,:,:)
1497	INTEGER, POINTER :: msk_glo2d(:,:)
1498	DOUBLE PRECISION, POINTER :: srf_glo2d(:,:)
1499	! Local variables
1500	INTEGER :: ierr
1501	CHARACTER(len=256) :: grid_filename
1502	CHARACTER(len=256) :: mask_filename
1503	CHARACTER(len=256) :: surf_filename
1504	CHARACTER(len=8) :: cl_nam ! cl_grd+.lon,+.lat ...
1505	INTEGER :: il_file_id, il_lon_id, il_lat_id, il_clo_id, il_cla_id
1506	INTEGER :: il_msk_id, il_srf_id
1507	INTEGER :: lon_dims, lat_dims
1508	INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: lon_dims_ids,lat_dims_ids
1509	INTEGER, DIMENSION(NF90_MAX_VAR_DIMS) :: lon_dims_len,lat_dims_len
1510	INTEGER :: i
1511	INTEGER, DIMENSION(:), ALLOCATABLE :: ila_dim, ila_what
1512
1513	grid_filename = TRIM(params%grids_file)
1514	mask_filename = TRIM(params%masks_file)
1515	surf_filename = TRIM(params%areas_file)
1516
1517	ierr = NF90_OPEN(TRIM(grid_filename), NF90_NOWRITE, il_file_id)
1518	cl_nam = TRIM(params%domain)//".lon"
1519	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_lon_id)
1520	ierr = NF90_INQUIRE_VARIABLE(il_file_id, varid=il_lon_id, ndims=lon_dims, dimids=lon_dims_ids)
1521	DO i=1,lon_dims
1522	ierr = NF90_INQUIRE_DIMENSION(ncid=il_file_id,dimid=lon_dims_ids(i),len=lon_dims_len(i))
1523	ENDDO
1524	ni_glo=lon_dims_len(1)
1525	nj_glo=lon_dims_len(2)
1526
1527	cl_nam = TRIM(params%domain)//".lat"
1528	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_lat_id)
1529	ierr = NF90_INQUIRE_VARIABLE(il_file_id, varid=il_lat_id, ndims=lat_dims, dimids=lat_dims_ids)
1530	DO i=1,lat_dims
1531	ierr = NF90_INQUIRE_DIMENSION(ncid=il_file_id,dimid=lat_dims_ids(i),len=lat_dims_len(i))
1532	ENDDO
1533	ALLOCATE(lon_glo2D(0:ni_glo-1,0:nj_glo-1), lat_glo2D(0:ni_glo-1,0:nj_glo-1))
1534	ALLOCATE(ila_what(2), ila_dim(2))
1535	ila_what(:)=1
1536	ila_dim(:)=[ni_glo, nj_glo]
1537	ierr = NF90_GET_VAR (il_file_id, il_lon_id, lon_glo2D(0:ni_glo-1,0:nj_glo-1), ila_what, ila_dim)
1538	ierr = NF90_GET_VAR (il_file_id, il_lat_id, lat_glo2D(0:ni_glo-1,0:nj_glo-1), ila_what, ila_dim)
1539	DEALLOCATE(ila_what, ila_dim)
1540
1541	cl_nam = TRIM(params%domain)//".clo"
1542	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_clo_id)
1543	ierr = NF90_INQUIRE_VARIABLE(il_file_id, varid=il_clo_id, ndims=lon_dims, dimids=lon_dims_ids)
1544	DO i=1,lon_dims
1545	ierr = NF90_INQUIRE_DIMENSION(ncid=il_file_id,dimid=lon_dims_ids(i),len=lon_dims_len(i))
1546	ENDDO
1547	ncrn=lon_dims_len(3)
1548
1549	cl_nam = TRIM(params%domain)//".cla"
1550	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_cla_id)
1551	ierr = NF90_INQUIRE_VARIABLE(il_file_id, varid=il_lat_id, ndims=lat_dims, dimids=lat_dims_ids)
1552	DO i=1,lat_dims
1553	ierr = NF90_INQUIRE_DIMENSION(ncid=il_file_id,dimid=lat_dims_ids(i),len=lat_dims_len(i))
1554	ENDDO
1555
1556	ALLOCATE(clo_glo2D(0:ni_glo-1,0:nj_glo-1,ncrn), cla_glo2D(0:ni_glo-1,0:nj_glo-1,ncrn))
1557	ALLOCATE(ila_what(3), ila_dim(3))
1558	ila_what(:)=1
1559	ila_dim(:)=[ni_glo, nj_glo, ncrn]
1560	ierr = NF90_GET_VAR (il_file_id, il_clo_id, clo_glo2D(0:ni_glo-1,0:nj_glo-1,1:ncrn), ila_what, ila_dim)
1561	ierr = NF90_GET_VAR (il_file_id, il_cla_id, cla_glo2D(0:ni_glo-1,0:nj_glo-1,1:ncrn), ila_what, ila_dim)
1562	DEALLOCATE(ila_what, ila_dim)
1563
1564	ierr = NF90_CLOSE(il_file_id)
1565
1566	ierr = NF90_OPEN(TRIM(mask_filename), NF90_NOWRITE, il_file_id)
1567	cl_nam = TRIM(params%domain)//".msk"
1568	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_msk_id)
1569	ALLOCATE(msk_glo2D(0:ni_glo-1,0:nj_glo-1))
1570	ALLOCATE(ila_what(2), ila_dim(2))
1571	ila_what(:)=1
1572	ila_dim(:)=[ni_glo, nj_glo]
1573	ierr = NF90_GET_VAR (il_file_id, il_msk_id, msk_glo2D(0:ni_glo-1,0:nj_glo-1), ila_what, ila_dim)
1574	DEALLOCATE(ila_what, ila_dim)
1575	ierr = NF90_CLOSE(il_file_id)
1576
1577	ierr = NF90_OPEN(TRIM(surf_filename), NF90_NOWRITE, il_file_id)
1578	cl_nam = TRIM(params%domain)//".srf"
1579	ierr = NF90_INQ_VARID(il_file_id, cl_nam, il_srf_id)
1580	ALLOCATE(srf_glo2D(0:ni_glo-1,0:nj_glo-1))
1581	ALLOCATE(ila_what(2), ila_dim(2))
1582	ila_what(:)=1
1583	ila_dim(:)=[ni_glo, nj_glo]
1584	ierr = NF90_GET_VAR (il_file_id, il_srf_id, srf_glo2D(0:ni_glo-1,0:nj_glo-1), ila_what, ila_dim)
1585	DEALLOCATE(ila_what, ila_dim)
1586	ierr = NF90_CLOSE(il_file_id)
1587
1588	END SUBROUTINE read_oasis_grid
1589
1590
1591
1592	SUBROUTINE get_decomposition(comm, params, domain_proc_size, domain_proc_rank, axis_proc_size, axis_proc_rank)
1593
1594	IMPLICIT NONE
1595	INTEGER,INTENT(IN) :: comm
1596	TYPE(tmodel_params) :: params
1597	INTEGER,INTENT(OUT) :: domain_proc_size
1598	INTEGER,INTENT(OUT) :: domain_proc_rank
1599	INTEGER,INTENT(OUT) :: axis_proc_size
1600	INTEGER,INTENT(OUT) :: axis_proc_rank
1601
1602	INTEGER :: mpi_rank,mpi_size,rank
1603	INTEGER :: ierr
1604	INTEGER :: n_domain,n_axis, min_dist, new_dist, best
1605	LOGICAL :: axis_layer
1606
1607	CALL MPI_COMM_RANK(comm,mpi_rank,ierr)
1608	CALL MPI_COMM_SIZE(comm,mpi_size,ierr)
1609
1610	IF (params%axis_proc_n > 0 ) THEN
1611	n_axis=params%axis_proc_n
1612	n_domain = mpi_size / n_axis
1613	axis_layer=.TRUE.
1614	ELSE IF (params%domain_proc_n > 0 ) THEN
1615	n_domain=params%domain_proc_n
1616	n_axis = mpi_size / n_domain
1617	axis_layer=.FALSE.
1618	ELSE
1619	IF (params%axis_proc_frac==0) THEN
1620	params%axis_proc_frac=1
1621	params%domain_proc_frac=mpi_size
1622	ELSE IF (params%domain_proc_frac==0) THEN
1623	params%domain_proc_frac=1
1624	params%axis_proc_frac=mpi_size
1625	ENDIF
1626
1627	n_domain = INT(SQRT(params%domain_proc_frac * mpi_size/ params%axis_proc_frac))
1628	n_axis = INT(SQRT(params%axis_proc_frac * mpi_size/ params%domain_proc_frac))
1629
1630
1631	min_dist= mpi_size - n_domain*n_axis
1632	best=0
1633
1634	new_dist = mpi_size -(n_domain+1)*n_axis
1635	IF (new_dist < min_dist .AND. new_dist >= 0 ) THEN
1636	min_dist=new_dist
1637	best=1
1638	ENDIF
1639
1640	new_dist=mpi_size-n_domain*(n_axis+1)
1641	IF (new_dist < min_dist .AND. new_dist >= 0 ) THEN
1642	min_dist=new_dist
1643	best=2
1644	ENDIF
1645
1646	IF (best==0) THEN
1647	ELSE IF (best==1) THEN
1648	n_domain=n_domain+1
1649	ELSE IF (best==2) THEN
1650	n_axis=n_axis+1
1651	ENDIF
1652
1653	IF ( MOD(mpi_size,n_axis) <= MOD(mpi_size,n_domain)) axis_layer=.TRUE.
1654
1655	ENDIF
1656
1657	IF ( axis_layer) THEN
1658	!!! n_axis layer
1659	IF (mpi_rank < MOD(mpi_size,n_axis)*(n_domain+1)) THEN
1660	axis_proc_rank=mpi_rank/(n_domain+1)
1661	domain_proc_rank=MOD(mpi_rank,n_domain+1)
1662	axis_proc_size=n_axis
1663	domain_proc_size=n_domain+1
1664	ELSE
1665	rank=mpi_rank-MOD(mpi_size,n_axis)*(n_domain+1)
1666	axis_proc_size=n_axis
1667	axis_proc_rank=MOD(mpi_size,n_axis)+rank/n_domain
1668	domain_proc_rank=MOD(rank,n_domain)
1669	domain_proc_size=n_domain
1670	ENDIF
1671	ELSE
1672	!!! n_domain column
1673	IF (mpi_rank < MOD(mpi_size,n_domain)*(n_axis+1)) THEN
1674	domain_proc_rank=mpi_rank/(n_axis+1)
1675	axis_proc_rank=MOD(mpi_rank,n_axis+1)
1676	domain_proc_size=n_domain
1677	axis_proc_size=n_axis+1
1678	ELSE
1679	rank=mpi_rank-MOD(mpi_size,n_domain)*(n_axis+1)
1680	domain_proc_size=n_domain
1681	domain_proc_rank=MOD(mpi_size,n_domain)+rank/n_axis
1682	axis_proc_rank=MOD(rank,n_axis)
1683	axis_proc_size=n_axis
1684	ENDIF
1685	ENDIF
1686
1687
1688	END SUBROUTINE get_decomposition
1689
1690	END PROGRAM oasis_testcase
1691

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

Context Navigation

source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/examples/regrid_environment/XIOS/oasis_testcase.f90 @ 6331

Download in other formats: