Context Navigation

source: codes/icosagcm/trunk/src/write_field.f90 @ 200

Last change on this file since 200 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM

File size: 62.8 KB

Line
1	module write_field
2	USE genmod
3	implicit none
4	PRIVATE
5	INTEGER,SAVE :: ncprec
6
7	TYPE ncvar
8	INTEGER :: size
9	INTEGER,POINTER :: nc_id(:)
10	INTEGER :: displ
11	END TYPE ncvar
12
13	INTEGER, PARAMETER :: MaxWriteField = 1000
14	INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId
15	TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId
16	INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex
17	CHARACTER(len=255), DIMENSION(MaxWriteField) :: FieldName
18
19	INTEGER,SAVE :: NbField = 0
20
21	PUBLIC init_writeField, writefield, close_files
22
23	CONTAINS
24
25	SUBROUTINE init_writeField
26	USE ioipsl
27	use netcdf_mod
28	IMPLICIT NONE
29	CHARACTER(LEN=255) :: netcdf_prec
30
31	netcdf_prec='float'
32	CALL getin("netcdf_prec",netcdf_prec)
33
34	SELECT CASE(TRIM(netcdf_prec))
35	CASE('float')
36	ncprec=NF90_FLOAT
37	CASE('double')
38	ncprec=NF90_DOUBLE
39	CASE default
40	PRINT*,'Bad selector for variable netcdf_prec : <', TRIM(netcdf_prec),"> options are <float>, <double>"
41	STOP
42	END SELECT
43	END SUBROUTINE init_writeField
44
45	function GetFieldIndex(name)
46	implicit none
47	integer :: GetFieldindex
48	character(len=*) :: name
49
50	character(len=255) :: TrueName
51	integer :: i
52
53
54	TrueName=TRIM(ADJUSTL(name))
55
56	GetFieldIndex=-1
57	do i=1,NbField
58	if (TrueName==FieldName(i)) then
59	GetFieldIndex=i
60	exit
61	endif
62	enddo
63	end function GetFieldIndex
64
65	SUBROUTINE Writefield(name_in,field,nind,once)
66	USE domain_mod
67	USE field_mod
68	USE transfert_mpi_mod
69	USE dimensions
70	USE mpipara
71	IMPLICIT NONE
72	CHARACTER(LEN=*),INTENT(IN) :: name_in
73	TYPE(t_field),POINTER :: field(:)
74	INTEGER,OPTIONAL,INTENT(IN) :: nind
75	LOGICAL,OPTIONAL,INTENT(IN) :: once
76	LOGICAL :: once_
77
78	TYPE(t_field),POINTER :: field_glo(:)
79
80	!$OMP MASTER
81	IF(PRESENT(once)) THEN
82	once_=once
83	ELSE
84	once_=.FALSE.
85	END IF
86
87	IF (field(1)%ndim==2) THEN
88	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type)
89	ELSE IF (field(1)%ndim==3) THEN
90	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3)
91	ELSE IF (field(1)%ndim==4) THEN
92	CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4)
93	ENDIF
94
95	CALL gather_field(field,field_glo)
96
97	IF (mpi_rank==0) THEN
98	IF (PRESENT(nind)) THEN
99	CALL writefield_gen(name_in,field_glo,domain_glo,nind,once=once_)
100	ELSE
101	CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo,once=once_)
102	ENDIF
103	ENDIF
104
105	CALL deallocate_field_glo(field_glo)
106	!$OMP END MASTER
107
108	END SUBROUTINE writefield
109
110	! SUBROUTINE Writefield(name_in,field,nind)
111	! USE netcdf
112	! USE domain_mod
113	! use field_mod
114	! USE dimensions
115	! USE geometry
116	! IMPLICIT NONE
117	! CHARACTER(LEN=*),INTENT(IN) :: name_in
118	! TYPE(t_field),POINTER :: field(:)
119	! INTEGER,OPTIONAL,INTENT(IN) :: nind
120	! REAL(r8),ALLOCATABLE :: field_val2d(:)
121	! REAL(r8),ALLOCATABLE :: field_val3d(:,:)
122	! REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
123	! TYPE(t_domain),POINTER :: d
124	! INTEGER :: Index
125	! INTEGER :: ind,i,j,k,n,ncell,q
126	! INTEGER :: iie,jje,iin,jjn
127	! INTEGER :: status
128	! CHARACTER(len=255) :: name
129	! CHARACTER(len=255) :: str_ind
130	! INTEGER :: ind_b,ind_e
131	! INTEGER :: halo_size
132	! LOGICAL :: single
133	!
134	!
135	! name=TRIM(ADJUSTL(name_in))
136
137	! IF (PRESENT(nind)) THEN
138	! name=TRIM(name)//"_"//TRIM(int2str(nind))
139	! PRINT *,"NAME",nind,int2str(nind),name
140	! ind_b=nind
141	! ind_e=nind
142	! halo_size=1
143	! single=.TRUE.
144	! ELSE
145	! ind_b=1
146	! ind_e=ndomain
147	! halo_size=0
148	! single=.FALSE.
149	! ENDIF
150
151	! Index=GetFieldIndex(name)
152	! if (Index==-1) then
153	! call create_header(name,field,nind)
154	! Index=GetFieldIndex(name)
155	! else
156	! FieldIndex(Index)=FieldIndex(Index)+1.
157	! endif
158	!
159	! IF (Field(ind_b)%field_type==field_T) THEN
160	! ncell=1
161	! DO ind=ind_b,ind_e
162	! d=>domain(ind)
163	! IF (Field(ind)%field_type/=field_T) THEN
164	! PRINT *,"Writefield, grille non geree"
165	! RETURN
166	! ENDIF
167
168	! n=0
169	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
170	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
171	! IF (d%own(i,j) .OR. single) n=n+1
172	! ENDDO
173	! ENDDO
174
175	! IF (field(ind)%ndim==2) THEN
176	! ALLOCATE(Field_val2d(n))
177	! n=0
178	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
179	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
180	! k=d%iim*(j-1)+i
181	! IF (d%own(i,j) .OR. single) THEN
182	! n=n+1
183	! Field_val2d(n)=field(ind)%rval2d(k)
184	! ENDIF
185	! ENDDO
186	! ENDDO
187	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
188	! start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
189	! DEALLOCATE(field_val2d)
190	! ELSE IF (field(ind)%ndim==3) THEN
191	! ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
192	! n=0
193	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
194	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
195	! k=d%iim*(j-1)+i
196	! IF (d%own(i,j) .OR. single) THEN
197	! n=n+1
198	! Field_val3d(n,:)=field(ind)%rval3d(k,:)
199	! ENDIF
200	! ENDDO
201	! ENDDO
202	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
203	! count=(/n,size(field(1)%rval3d,2),1 /))
204	! DEALLOCATE(field_val3d)
205	! ELSE IF (field(1)%ndim==4) THEN
206
207	! DO q=1,FieldVarId(index)%size
208	!
209	! ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
210	! n=0
211	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
212	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
213	! k=d%iim*(j-1)+i
214	! IF (d%own(i,j) .OR. single) THEN
215	! n=n+1
216	! Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
217	! ENDIF
218	! ENDDO
219	! ENDDO
220	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
221	! count=(/n,size(field(1)%rval4d,2),1 /))
222	! DEALLOCATE(field_val3d)
223	! ENDDO
224	! ENDIF
225	!
226	! PRINT *,NF90_STRERROR(status)
227	! ncell=ncell+n
228
229	! ENDDO
230	!
231	! ELSE IF (Field(ind_b)%field_type==field_Z) THEN
232	! ncell=1
233	! n=0
234	! DO ind=ind_b,ind_e
235	! d=>domain(ind)
236	! CALL swap_geometry(ind)
237	! CALL swap_dimensions(ind)
238	!
239	! n=0
240	! DO j=jj_begin+1,jj_end
241	! DO i=ii_begin,ii_end-1
242	! n=n+1
243	! ENDDO
244	! ENDDO
245
246	! DO j=jj_begin,jj_end-1
247	! DO i=ii_begin+1,ii_end
248	! n=n+1
249	! ENDDO
250	! ENDDO
251
252	! IF (field(ind)%ndim==2) THEN
253	! ALLOCATE(Field_val2d(n))
254
255	! n=0
256	! DO j=jj_begin+1,jj_end
257	! DO i=ii_begin,ii_end-1
258	! n=n+1
259	! k=iim*(j-1)+i
260	! Field_val2d(n)=field(ind)%rval2d(k+z_down)
261	! ENDDO
262	! ENDDO
263
264	! DO j=jj_begin,jj_end-1
265	! DO i=ii_begin+1,ii_end
266	! n=n+1
267	! k=iim*(j-1)+i
268	! Field_val2d(n)=field(ind)%rval2d(k+z_up)
269	! ENDDO
270	! ENDDO
271
272	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
273	! Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
274	! DEALLOCATE(field_val2d)
275
276	! ELSE IF (field(ind)%ndim==3) THEN
277	! ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
278	! n=0
279	! DO j=jj_begin+1,jj_end
280	! DO i=ii_begin,ii_end-1
281	! n=n+1
282	! k=iim*(j-1)+i
283	! Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
284	! ENDDO
285	! ENDDO
286
287	! DO j=jj_begin,jj_end-1
288	! DO i=ii_begin+1,ii_end
289	! n=n+1
290	! k=iim*(j-1)+i
291	! Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
292	! ENDDO
293	! ENDDO
294	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /), &
295	! count=(/n,size(field(1)%rval3d,2),1 /))
296	! DEALLOCATE(field_val3d)
297	! ELSE IF (field(1)%ndim==4) THEN
298
299	! DO q=1,FieldVarId(index)%size
300	! ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
301	! n=0
302	! DO j=jj_begin+1,jj_end
303	! DO i=ii_begin,ii_end-1
304	! n=n+1
305	! k=iim*(j-1)+i
306	! Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
307	! ENDDO
308	! ENDDO
309
310	! DO j=jj_begin,jj_end-1
311	! DO i=ii_begin+1,ii_end
312	! n=n+1
313	! k=iim*(j-1)+i
314	! Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
315	! ENDDO
316	! ENDDO
317
318	! status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /), &
319	! count=(/n,size(field(1)%rval4d,2),1 /))
320	! DEALLOCATE(field_val3d)
321	! ENDDO
322	! ENDIF
323	!
324	! PRINT *,NF90_STRERROR(status)
325	! ncell=ncell+n
326
327	! ENDDO
328	!
329	! ENDIF
330	! status=NF90_SYNC(FieldId(Index))
331	!
332	! END SUBROUTINE Writefield
333
334
335	SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once )
336	USE netcdf_mod
337	USE domain_mod
338	USE field_mod
339	IMPLICIT NONE
340	CHARACTER(LEN=*),INTENT(IN) :: name_in
341	TYPE(t_field), POINTER :: field(:)
342	TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
343	INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
344	INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
345	REAL(r8),ALLOCATABLE :: field_val2d(:)
346	REAL(r8),ALLOCATABLE :: field_val3d(:,:)
347	REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
348	LOGICAL, INTENT(IN) :: once
349	TYPE(t_domain),POINTER :: d
350	INTEGER :: Index
351	INTEGER :: ind,i,j,k,n,ncell,q
352	INTEGER :: iie,jje,iin,jjn
353	INTEGER :: status
354	CHARACTER(len=255) :: name
355	CHARACTER(len=255) :: str_ind
356	INTEGER :: ind_b,ind_e
357	INTEGER :: halo_size
358	LOGICAL :: single
359
360	name=TRIM(ADJUSTL(name_in))
361
362	IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
363	name=TRIM(name)//"_"//TRIM(int2str(ind_b))
364	ind_b=ind_b_in
365	ind_e=ind_b_in
366	halo_size=1
367	single=.TRUE.
368	ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
369	name=TRIM(name)//"_"//TRIM(int2str(ind_e))
370	ind_b=ind_e_in
371	ind_e=ind_e_in
372	halo_size=1
373	single=.TRUE.
374	ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
375	ind_b=ind_b_in
376	ind_e=ind_e_in
377	halo_size=0
378	single=.FALSE.
379	ELSE
380	halo_size=0
381	ind_b=1
382	ind_e=ndomain
383	single=.FALSE.
384	ENDIF
385
386	Index=GetFieldIndex(name)
387	if (Index==-1) then
388	call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
389	Index=GetFieldIndex(name)
390	else
391	FieldIndex(Index)=FieldIndex(Index)+1.
392	endif
393
394	IF (Field(ind_b)%field_type==field_T) THEN
395
396	ncell=0
397	DO ind=ind_b,ind_e
398	d=>domain_type(ind)
399	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
400	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
401	IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1
402	ENDDO
403	ENDDO
404	ENDDO
405
406	IF (field(1)%ndim==2) THEN
407	ALLOCATE(Field_val2d(ncell))
408	n=0
409	DO ind=ind_b,ind_e
410	d=>domain_type(ind)
411	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
412	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
413	k=d%iim*(j-1)+i
414	IF (d%assign_domain(i,j)==ind .OR. single) THEN
415	n=n+1
416	Field_val2d(n)=field(ind)%rval2d(k)
417	ENDIF
418	ENDDO
419	ENDDO
420	ENDDO
421
422	IF (once) THEN
423	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
424	start=(/ 1 /),count=(/ncell /))
425	ELSE
426	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
427	start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
428	ENDIF
429
430	DEALLOCATE(field_val2d)
431	ELSE IF (field(1)%ndim==3) THEN
432	ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
433	n=0
434	DO ind=ind_b,ind_e
435	d=>domain_type(ind)
436	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
437	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
438	k=d%iim*(j-1)+i
439	IF (d%assign_domain(i,j)==ind .OR. single) THEN
440	n=n+1
441	Field_val3d(n,:)=field(ind)%rval3d(k,:)
442	ENDIF
443	ENDDO
444	ENDDO
445	ENDDO
446
447	PRINT *, 'Writefield ', TRIM(name), MAXVAL(Field_val3d(:,1)), MINVAL(Field_val3d(:,1)) ! FIXME
448
449	IF (once) THEN
450	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /), &
451	count=(/ncell,size(field(1)%rval3d,2) /))
452	ELSE
453	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
454	count=(/ncell,size(field(1)%rval3d,2),1 /))
455	ENDIF
456
457	DEALLOCATE(field_val3d)
458	ELSE IF (field(1)%ndim==4) THEN
459
460	DO q=1,FieldVarId(index)%size
461
462	ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
463	n=0
464	DO ind=ind_b,ind_e
465	d=>domain_type(ind)
466	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
467	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
468	k=d%iim*(j-1)+i
469	IF (d%assign_domain(i,j)==ind .OR. single) THEN
470	n=n+1
471	Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
472	ENDIF
473	ENDDO
474	ENDDO
475	ENDDO
476
477	IF (once) THEN
478	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1 /), &
479	count=(/ncell,size(field(1)%rval4d,2),1 /))
480	ELSE
481	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
482	count=(/ncell,size(field(1)%rval4d,2),1 /))
483	ENDIF
484	DEALLOCATE(field_val3d)
485	ENDDO
486	ENDIF
487
488
489	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
490
491	ncell=0
492	DO ind=ind_b,ind_e
493	d=>domain_type(ind)
494
495	DO j=d%jj_begin+1,d%jj_end
496	DO i=d%ii_begin,d%ii_end-1
497	ncell=ncell+1
498	ENDDO
499	ENDDO
500
501	DO j=d%jj_begin,d%jj_end-1
502	DO i=d%ii_begin+1,d%ii_end
503	ncell=ncell+1
504	ENDDO
505	ENDDO
506	ENDDO
507
508	IF (field(1)%ndim==2) THEN
509	ALLOCATE(Field_val2d(ncell))
510
511	n=0
512
513	DO ind=ind_b,ind_e
514	d=>domain_type(ind)
515	DO j=d%jj_begin+1,d%jj_end
516	DO i=d%ii_begin,d%ii_end-1
517	n=n+1
518	k=d%iim*(j-1)+i
519	Field_val2d(n)=field(ind)%rval2d(k+d%z_down)
520	ENDDO
521	ENDDO
522
523	DO j=d%jj_begin,d%jj_end-1
524	DO i=d%ii_begin+1,d%ii_end
525	n=n+1
526	k=d%iim*(j-1)+i
527	Field_val2d(n)=field(ind)%rval2d(k+d%z_up)
528	ENDDO
529	ENDDO
530	ENDDO
531
532	IF (once) THEN
533	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
534	Field_val2d,start=(/ 1 /),count=(/ncell /))
535	ELSE
536	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
537	Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
538	ENDIF
539	DEALLOCATE(field_val2d)
540
541	ELSE IF (field(1)%ndim==3) THEN
542	ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
543	n=0
544	DO ind=ind_b,ind_e
545	d=>domain_type(ind)
546	DO j=d%jj_begin+1,d%jj_end
547	DO i=d%ii_begin,d%ii_end-1
548	n=n+1
549	k=d%iim*(j-1)+i
550	Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:)
551	ENDDO
552	ENDDO
553
554	DO j=d%jj_begin,d%jj_end-1
555	DO i=d%ii_begin+1,d%ii_end
556	n=n+1
557	k=d%iim*(j-1)+i
558	Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:)
559	ENDDO
560	ENDDO
561	ENDDO
562
563	IF (once) THEN
564	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /), &
565	count=(/ncell,size(field(1)%rval3d,2) /))
566	ELSE
567	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /), &
568	count=(/ncell,size(field(1)%rval3d,2),1 /))
569	ENDIF
570	DEALLOCATE(field_val3d)
571
572	ELSE IF (field(1)%ndim==4) THEN
573
574	DO q=1,FieldVarId(index)%size
575	ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
576	n=0
577	DO ind=ind_b,ind_e
578	d=>domain_type(ind)
579	DO j=d%jj_begin+1,d%jj_end
580	DO i=d%ii_begin,d%ii_end-1
581	n=n+1
582	k=d%iim*(j-1)+i
583	Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q)
584	ENDDO
585	ENDDO
586
587	DO j=d%jj_begin,d%jj_end-1
588	DO i=d%ii_begin+1,d%ii_end
589	n=n+1
590	k=d%iim*(j-1)+i
591	Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q)
592	ENDDO
593	ENDDO
594	ENDDO
595
596	IF (once) THEN
597	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1 /), &
598	count=(/ncell,size(field(1)%rval4d,2) /))
599	ELSE
600	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /), &
601	count=(/ncell,size(field(1)%rval4d,2),1 /))
602	ENDIF
603	DEALLOCATE(field_val3d)
604	ENDDO
605	ENDIF
606
607	ENDIF
608	status=NF90_SYNC(FieldId(Index))
609
610	END SUBROUTINE Writefield_gen
611
612
613
614	SUBROUTINE Writefield_mpi(name_in,field,nind)
615	USE netcdf_mod
616	USE domain_mod
617	use field_mod
618	USE dimensions
619	USE geometry
620	IMPLICIT NONE
621	CHARACTER(LEN=*),INTENT(IN) :: name_in
622	TYPE(t_field),POINTER :: field(:)
623	INTEGER,OPTIONAL,INTENT(IN) :: nind
624	REAL(r8),ALLOCATABLE :: field_val2d(:)
625	REAL(r8),ALLOCATABLE :: field_val3d(:,:)
626	REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
627	TYPE(t_domain),POINTER :: d
628	INTEGER :: Index
629	INTEGER :: ind,i,j,l,k,n,ncell,q
630	INTEGER :: iie,jje,iin,jjn
631	INTEGER :: status
632	CHARACTER(len=255) :: name
633	CHARACTER(len=255) :: str_ind
634	INTEGER :: ind_b,ind_e
635	INTEGER :: halo_size
636	LOGICAL :: single
637	INTEGER :: displ
638
639
640	name=TRIM(ADJUSTL(name_in))
641
642	IF (PRESENT(nind)) THEN
643	name=TRIM(name)//"_"//TRIM(int2str(nind))
644	PRINT *,"NAME",nind,int2str(nind),name
645	ind_b=nind
646	ind_e=nind
647	halo_size=1
648	single=.TRUE.
649	ELSE
650	ind_b=1
651	ind_e=ndomain
652	halo_size=0
653	single=.FALSE.
654	ENDIF
655
656	Index=GetFieldIndex(name)
657	if (Index==-1) then
658	call create_header_mpi(name,field,nind)
659	Index=GetFieldIndex(name)
660	else
661	FieldIndex(Index)=FieldIndex(Index)+1.
662	endif
663
664	IF (Field(ind_b)%field_type==field_T) THEN
665	ncell=1
666	DO ind=ind_b,ind_e
667	d=>domain(ind)
668	IF (Field(ind)%field_type/=field_T) THEN
669	PRINT *,"Writefield, grille non geree"
670	RETURN
671	ENDIF
672
673	n=0
674	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
675	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
676	IF (d%own(i,j) .OR. single) n=n+1
677	ENDDO
678	ENDDO
679
680	displ=FieldVarId(index)%displ
681
682	IF (field(ind)%ndim==2) THEN
683	ALLOCATE(Field_val2d(n))
684	n=0
685	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
686	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
687	k=d%iim*(j-1)+i
688	IF (d%own(i,j) .OR. single) THEN
689	n=n+1
690	Field_val2d(n)=field(ind)%rval2d(k)
691	ENDIF
692	ENDDO
693	ENDDO
694	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d, &
695	start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
696	DEALLOCATE(field_val2d)
697	ELSE IF (field(ind)%ndim==3) THEN
698	ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
699	n=0
700	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
701	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
702	k=d%iim*(j-1)+i
703	IF (d%own(i,j) .OR. single) THEN
704	n=n+1
705	Field_val3d(n,:)=field(ind)%rval3d(k,:)
706	ENDIF
707	ENDDO
708	ENDDO
709	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d, &
710	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
711	DEALLOCATE(field_val3d)
712
713	ELSE IF (field(1)%ndim==4) THEN
714
715	DO q=1,FieldVarId(index)%size
716
717	ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
718	n=0
719	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
720	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
721	k=d%iim*(j-1)+i
722	IF (d%own(i,j) .OR. single) THEN
723	n=n+1
724	Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
725	ENDIF
726	ENDDO
727	ENDDO
728
729	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
730	start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
731	DEALLOCATE(field_val3d)
732	ENDDO
733	ENDIF
734
735	ncell=ncell+n
736
737	ENDDO
738
739	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
740	ncell=1
741	n=0
742	DO ind=ind_b,ind_e
743	d=>domain(ind)
744	CALL swap_geometry(ind)
745	CALL swap_dimensions(ind)
746
747	n=0
748	DO j=jj_begin+1,jj_end
749	DO i=ii_begin,ii_end-1
750	n=n+1
751	ENDDO
752	ENDDO
753
754	DO j=jj_begin,jj_end-1
755	DO i=ii_begin+1,ii_end
756	n=n+1
757	ENDDO
758	ENDDO
759
760	displ=FieldVarId(index)%displ
761
762	IF (field(ind)%ndim==2) THEN
763	ALLOCATE(Field_val2d(n))
764
765	n=0
766	DO j=jj_begin+1,jj_end
767	DO i=ii_begin,ii_end-1
768	n=n+1
769	k=iim*(j-1)+i
770	Field_val2d(n)=field(ind)%rval2d(k+z_down)
771	ENDDO
772	ENDDO
773
774	DO j=jj_begin,jj_end-1
775	DO i=ii_begin+1,ii_end
776	n=n+1
777	k=iim*(j-1)+i
778	Field_val2d(n)=field(ind)%rval2d(k+z_up)
779	ENDDO
780	ENDDO
781
782	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
783	Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
784	DEALLOCATE(field_val2d)
785
786	ELSE IF (field(ind)%ndim==3) THEN
787	ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
788	n=0
789	DO j=jj_begin+1,jj_end
790	DO i=ii_begin,ii_end-1
791	n=n+1
792	k=iim*(j-1)+i
793	Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
794	ENDDO
795	ENDDO
796
797	DO j=jj_begin,jj_end-1
798	DO i=ii_begin+1,ii_end
799	n=n+1
800	k=iim*(j-1)+i
801	Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
802	ENDDO
803	ENDDO
804
805	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d, &
806	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
807	DEALLOCATE(field_val3d)
808
809	ELSE IF (field(1)%ndim==4) THEN
810
811	DO q=1,FieldVarId(index)%size
812	ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
813	n=0
814	DO j=jj_begin+1,jj_end
815	DO i=ii_begin,ii_end-1
816	n=n+1
817	k=iim*(j-1)+i
818	Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
819	ENDDO
820	ENDDO
821
822	DO j=jj_begin,jj_end-1
823	DO i=ii_begin+1,ii_end
824	n=n+1
825	k=iim*(j-1)+i
826	Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
827	ENDDO
828	ENDDO
829
830	status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d, &
831	start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
832	DEALLOCATE(field_val3d)
833	ENDDO
834	ENDIF
835
836	ncell=ncell+n
837
838	ENDDO
839
840	ENDIF
841	status=NF90_SYNC(FieldId(Index))
842
843	END SUBROUTINE Writefield_mpi
844
845
846
847	! SUBROUTINE Create_header(name,field,nind)
848	! USE netcdf
849	! USE field_mod
850	! USE domain_mod
851	! USE spherical_geom_mod
852	! USE dimensions
853	! USE geometry
854	! IMPLICIT NONE
855	! CHARACTER(LEN=*) :: name
856	! TYPE(t_field),POINTER :: field(:)
857	! INTEGER,OPTIONAL,INTENT(IN) :: nind
858	! INTEGER :: ncell
859	! INTEGER :: nvert
860	! REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
861	! TYPE(t_domain),POINTER :: d
862	! INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
863	! INTEGER :: dim3id,dim4id
864	! INTEGER :: status
865	! INTEGER :: ind,i,j,k,n,q
866	! INTEGER :: iie,jje,iin,jjn
867	! INTEGER :: ind_b,ind_e
868	! INTEGER :: halo_size
869	! LOGICAL :: single
870	! INTEGER :: nij
871	!
872	! NbField=NbField+1
873	! FieldName(NbField)=TRIM(ADJUSTL(name))
874	! FieldIndex(NbField)=1
875	!
876	! IF (PRESENT(nind)) THEN
877	! ind_b=nind
878	! ind_e=nind
879	! halo_size=1
880	! single=.TRUE.
881	! ELSE
882	! ind_b=1
883	! ind_e=ndomain
884	! halo_size=0
885	! single=.FALSE.
886	! ENDIF
887	!
888	! ncell=0
889	!
890	! IF (Field(ind_b)%field_type==field_T) THEN
891	! nvert=6
892	!
893	! DO ind=ind_b,ind_e
894	! d=>domain(ind)
895	!
896	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
897	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
898	! IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
899	! ENDDO
900	! ENDDO
901
902	! END DO
903	!
904	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
905	! FieldId(NbField)=ncid
906	! status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
907	! status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
908
909	! IF (Field(ind_b)%ndim==2) THEN
910	! FieldVarId(NbField)%size=1
911	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
912	! ELSE IF (Field(ind_b)%ndim==3) THEN
913	! FieldVarId(NbField)%size=1
914	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
915	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
916	! ELSE IF (Field(1)%ndim==4) THEN
917	! FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
918	! ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
919	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
920	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
921	! ENDIF
922	!
923	! status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
924	!
925	! status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
926	! status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
927	! status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
928	! status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
929	! status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
930	! status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
931	! status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
932	! status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
933	! status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
934	! status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
935
936	! IF (Field(ind_b)%ndim==2) THEN
937	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
938	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
939	! ELSE IF (Field(ind_b)%ndim==3) THEN
940	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
941	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
942	! ELSE IF (Field(ind_b)%ndim==4) THEN
943	! DO i=1,FieldVarId(NbField)%size
944	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /), &
945	! FieldVarId(NbField)%nc_id(i))
946	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
947	! ENDDO
948	! ENDIF
949	!
950	!
951	! status = NF90_ENDDEF(ncid)
952
953	! ncell=1
954	! DO ind=ind_b,ind_e
955	! d=>domain(ind)
956	!
957	! n=0
958	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
959	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
960	! IF (single .OR. d%own(i,j)) n=n+1
961	! ENDDO
962	! ENDDO
963	!
964	! ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
965	!
966	! n=0
967	! DO j=d%jj_begin-halo_size,d%jj_end+halo_size
968	! DO i=d%ii_begin-halo_size,d%ii_end+halo_size
969	! IF (d%own(i,j) .OR. single) THEN
970	! n=n+1
971	! CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
972	! lon(n)=lon(n)*180/Pi
973	! lat(n)=lat(n)*180/Pi
974	! DO k=0,5
975	! CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
976	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
977	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
978	! ENDDO
979	! ENDIF
980	! ENDDO
981	! ENDDO
982	! status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
983	! status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
984	! status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
985	! status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
986	!
987	! ncell=ncell+n
988	! DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
989	! END DO
990
991	! ELSE IF (Field(ind_b)%field_type==field_Z) THEN
992	! nvert=3
993	! DO ind=ind_b,ind_e
994	! d=>domain(ind)
995	!
996	! DO j=d%jj_begin+1,d%jj_end
997	! DO i=d%ii_begin,d%ii_end-1
998	! ncell=ncell+1
999	! ENDDO
1000	! ENDDO
1001
1002	! DO j=d%jj_begin,d%jj_end-1
1003	! DO i=d%ii_begin+1,d%ii_end
1004	! ncell=ncell+1
1005	! ENDDO
1006	! ENDDO
1007
1008	! END DO
1009	!
1010	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1011	! FieldId(NbField)=ncid
1012	! status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
1013	! status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1014
1015	! IF (Field(ind_b)%ndim==2) THEN
1016	! FieldVarId(NbField)%size=1
1017	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
1018	! ELSE IF (Field(ind_b)%ndim==3) THEN
1019	! FieldVarId(NbField)%size=1
1020	! ALLOCATE(FieldVarId(NbField)%nc_id(1))
1021	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
1022	! ELSE IF (Field(1)%ndim==4) THEN
1023	! FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1024	! ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1025	! status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
1026	! ENDIF
1027
1028
1029	!
1030	! status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1031	!
1032	! status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1033	! status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1034	! status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1035	! status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1036	! status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1037	! status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1038	! status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1039	! status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1040	! status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1041	! status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1042
1043
1044	! IF (Field(ind_b)%ndim==2) THEN
1045	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1046	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1047	! ELSE IF (Field(ind_b)%ndim==3) THEN
1048	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1049	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1050	! ELSE IF (Field(ind_b)%ndim==4) THEN
1051	! DO q=1,FieldVarId(NbField)%size
1052	! status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE, &
1053	! (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1054	! status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
1055	! ENDDO
1056	! ENDIF
1057	!
1058	! status = NF90_ENDDEF(ncid)
1059
1060	! ncell=1
1061	! DO ind=ind_b,ind_e
1062	! d=>domain(ind)
1063	! CALL swap_geometry(ind)
1064	! CALL swap_dimensions(ind)
1065	!
1066	! n=0
1067	! DO j=jj_begin+1,jj_end
1068	! DO i=ii_begin,ii_end-1
1069	! n=n+1
1070	! ENDDO
1071	! ENDDO
1072
1073	! DO j=jj_begin,jj_end-1
1074	! DO i=ii_begin+1,ii_end
1075	! n=n+1
1076	! ENDDO
1077	! ENDDO
1078
1079	! ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1080	!
1081	! n=0
1082	!
1083	! DO j=jj_begin+1,jj_end
1084	! DO i=ii_begin,ii_end-1
1085	! nij=(j-1)*iim+i
1086	! n=n+1
1087	! CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
1088	! lon(n)=lon(n)*180/Pi
1089	!! IF (lon(n)<0) lon(n)=lon(n)+360
1090	! lat(n)=lat(n)*180/Pi
1091	! CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1092	! CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1093	! CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1094
1095	! DO k=0,2
1096	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1097	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1098	! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
1099	! ENDDO
1100	! ENDDO
1101	! ENDDO
1102
1103	! DO j=jj_begin,jj_end-1
1104	! DO i=ii_begin+1,ii_end
1105	! nij=(j-1)*iim+i
1106	! n=n+1
1107	! CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
1108	! lon(n)=lon(n)*180/Pi
1109	! IF (lon(n)<0) lon(n)=lon(n)+360
1110	! lat(n)=lat(n)*180/Pi
1111	! CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1112	! CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1113	! CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1114
1115	! DO k=0,2
1116	! bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1117	! bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1118	! IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
1119	! ENDDO
1120	! ENDDO
1121	! ENDDO
1122	!
1123	!
1124	! status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
1125	! status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
1126	! status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
1127	! status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
1128	!
1129	! ncell=ncell+n
1130	! DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1131	! END DO
1132	! ENDIF
1133
1134
1135	!
1136	! status = NF90_CLOSE(ncid)
1137
1138	! END SUBROUTINE Create_Header
1139
1140
1141
1142
1143	SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
1144	USE netcdf_mod
1145	USE field_mod
1146	USE domain_mod
1147	USE metric
1148	USE spherical_geom_mod
1149	IMPLICIT NONE
1150	CHARACTER(LEN=*),INTENT(IN) :: name_in
1151	TYPE(t_field),POINTER :: field(:)
1152	TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
1153	INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
1154	INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
1155	LOGICAL,INTENT(IN) :: once
1156
1157	INTEGER :: ncell
1158	INTEGER :: nvert
1159	REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
1160	TYPE(t_domain),POINTER :: d
1161	INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
1162	INTEGER :: dim3id,dim4id
1163	INTEGER :: status
1164	INTEGER :: ind,i,j,k,n,q
1165	INTEGER :: iie,jje,iin,jjn
1166	INTEGER :: ind_b,ind_e
1167	INTEGER :: halo_size
1168	LOGICAL :: single
1169	INTEGER :: nij
1170	CHARACTER(LEN=255) :: name
1171
1172
1173	name=TRIM(ADJUSTL(name_in))
1174
1175	IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
1176	name=TRIM(name)//"_"//TRIM(int2str(ind_b))
1177	ind_b=ind_b_in
1178	ind_e=ind_b_in
1179	halo_size=1
1180	single=.TRUE.
1181	ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
1182	name=TRIM(name)//"_"//TRIM(int2str(ind_e))
1183	ind_b=ind_e_in
1184	ind_e=ind_e_in
1185	halo_size=1
1186	single=.TRUE.
1187	ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
1188	ind_b=ind_b_in
1189	ind_e=ind_e_in
1190	halo_size=0
1191	single=.FALSE.
1192	ELSE
1193	halo_size=0
1194	ind_b=1
1195	ind_e=ndomain
1196	single=.FALSE.
1197	ENDIF
1198
1199	NbField=NbField+1
1200	FieldName(NbField)=TRIM(ADJUSTL(name))
1201	FieldIndex(NbField)=1
1202
1203	ncell=0
1204
1205	IF (Field(ind_b)%field_type==field_T) THEN
1206	nvert=6
1207
1208	DO ind=ind_b,ind_e
1209	d=>domain_type(ind)
1210
1211	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1212	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1213	IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1
1214	ENDDO
1215	ENDDO
1216
1217	END DO
1218
1219	status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1220	FieldId(NbField)=ncid
1221	status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId)
1222	status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid)
1223
1224	IF (Field(ind_b)%ndim==2) THEN
1225	FieldVarId(NbField)%size=1
1226	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1227	ELSE IF (Field(ind_b)%ndim==3) THEN
1228	FieldVarId(NbField)%size=1
1229	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1230	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1231	ELSE IF (Field(1)%ndim==4) THEN
1232	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1233	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1234	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1235	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
1236	ENDIF
1237
1238	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1239
1240	status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId)
1241	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1242	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1243	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i")
1244	status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
1245	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1246	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1247	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i")
1248	status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1249	status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1250
1251	IF (Field(ind_b)%ndim==2) THEN
1252	IF (once) THEN
1253	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
1254	ELSE
1255	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1256	ENDIF
1257	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
1258	ELSE IF (Field(ind_b)%ndim==3) THEN
1259	IF (once) THEN
1260	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
1261	ELSE
1262	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1263	ENDIF
1264	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
1265	ELSE IF (Field(ind_b)%ndim==4) THEN
1266	DO i=1,FieldVarId(NbField)%size
1267	IF (once) THEN
1268	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /), &
1269	FieldVarId(NbField)%nc_id(i))
1270	ELSE
1271	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /), &
1272	FieldVarId(NbField)%nc_id(i))
1273	ENDIF
1274	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
1275	ENDDO
1276	ENDIF
1277
1278
1279	status = NF90_ENDDEF(ncid)
1280
1281	ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
1282
1283	n=0
1284	DO ind=ind_b,ind_e
1285	d=>domain_type(ind)
1286	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1287	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1288	IF (d%assign_domain(i,j)==ind .OR. single) THEN
1289	n=n+1
1290	CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
1291	lon(n)=lon(n)*180/Pi
1292	lat(n)=lat(n)*180/Pi
1293	DO k=0,5
1294	CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
1295	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1296	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1297	ENDDO
1298	ENDIF
1299	ENDDO
1300	ENDDO
1301	ENDDO
1302	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
1303	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
1304	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1305	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1306
1307	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1308
1309	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
1310	nvert=3
1311	DO ind=ind_b,ind_e
1312	d=>domain_type(ind)
1313
1314	DO j=d%jj_begin+1,d%jj_end
1315	DO i=d%ii_begin,d%ii_end-1
1316	ncell=ncell+1
1317	ENDDO
1318	ENDDO
1319
1320	DO j=d%jj_begin,d%jj_end-1
1321	DO i=d%ii_begin+1,d%ii_end
1322	ncell=ncell+1
1323	ENDDO
1324	ENDDO
1325
1326	END DO
1327
1328	status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1329	FieldId(NbField)=ncid
1330	status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId)
1331	status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid)
1332
1333	IF (Field(ind_b)%ndim==2) THEN
1334	FieldVarId(NbField)%size=1
1335	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1336	ELSE IF (Field(ind_b)%ndim==3) THEN
1337	FieldVarId(NbField)%size=1
1338	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1339	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1340	ELSE IF (Field(1)%ndim==4) THEN
1341	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1342	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1343	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1344	ENDIF
1345
1346
1347
1348	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1349
1350	status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId)
1351	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1352	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1353	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v")
1354	status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
1355	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1356	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1357	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v")
1358	status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1359	status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1360
1361
1362	IF (Field(ind_b)%ndim==2) THEN
1363	IF (once) THEN
1364	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
1365	ELSE
1366	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1367	ENDIF
1368	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
1369	ELSE IF (Field(ind_b)%ndim==3) THEN
1370	IF (once) THEN
1371	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
1372	ELSE
1373	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1374	ENDIF
1375	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
1376	ELSE IF (Field(ind_b)%ndim==4) THEN
1377	DO q=1,FieldVarId(NbField)%size
1378	IF (once) THEN
1379	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec, &
1380	(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(q))
1381	ELSE
1382	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec, &
1383	(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1384	ENDIF
1385	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon_v lat_v")
1386	ENDDO
1387	ENDIF
1388
1389	status = NF90_ENDDEF(ncid)
1390
1391	ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
1392
1393	n=0
1394
1395	DO ind=ind_b,ind_e
1396	d=>domain_type(ind)
1397	DO j=d%jj_begin+1,d%jj_end
1398	DO i=d%ii_begin,d%ii_end-1
1399	nij=(j-1)*d%iim+i
1400	n=n+1
1401	CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n))
1402	lon(n)=lon(n)*180/Pi
1403	lat(n)=lat(n)*180/Pi
1404
1405	CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
1406	CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n))
1407	CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n))
1408
1409	DO k=0,2
1410	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1411	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1412	ENDDO
1413	ENDDO
1414	ENDDO
1415
1416	DO j=d%jj_begin,d%jj_end-1
1417	DO i=d%ii_begin+1,d%ii_end
1418	nij=(j-1)*d%iim+i
1419	n=n+1
1420	CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n))
1421	lon(n)=lon(n)*180/Pi
1422	lat(n)=lat(n)*180/Pi
1423	CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
1424	CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n))
1425	CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n))
1426
1427	DO k=0,2
1428	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1429	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1430	ENDDO
1431	ENDDO
1432	ENDDO
1433	ENDDO
1434
1435	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
1436	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
1437	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1438	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
1439
1440	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1441	ENDIF
1442
1443
1444	END SUBROUTINE Create_Header_gen
1445
1446	SUBROUTINE Create_header_mpi(name,field,nind)
1447	USE netcdf_mod
1448	USE field_mod
1449	USE domain_mod
1450	USE spherical_geom_mod
1451	USE dimensions
1452	USE geometry
1453	USE mpi_mod
1454	USE mpipara
1455	IMPLICIT NONE
1456	CHARACTER(LEN=*) :: name
1457	CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj
1458	TYPE(t_field),POINTER :: field(:)
1459	INTEGER,OPTIONAL,INTENT(IN) :: nind
1460	INTEGER :: ncell
1461	INTEGER :: nvert
1462	REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
1463	TYPE(t_domain),POINTER :: d
1464	INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
1465	INTEGER :: dim3id,dim4id
1466	INTEGER :: status
1467	INTEGER :: ind,i,j,k,n,q
1468	INTEGER :: iie,jje,iin,jjn
1469	INTEGER :: ind_b,ind_e
1470	INTEGER :: halo_size
1471	LOGICAL :: single
1472	INTEGER :: nij
1473	INTEGER :: ncell_glo(0:mpi_size-1)
1474	INTEGER :: displ, ncell_tot
1475
1476
1477	NbField=NbField+1
1478	name_adj=TRIM(ADJUSTL(name)) ! work around ICE with pgf90
1479	FieldName(NbField)=name_adj
1480	FieldIndex(NbField)=1
1481
1482	IF (PRESENT(nind)) THEN
1483	ind_b=nind
1484	ind_e=nind
1485	halo_size=1
1486	single=.TRUE.
1487	ELSE
1488	ind_b=1
1489	ind_e=ndomain
1490	halo_size=0
1491	single=.FALSE.
1492	ENDIF
1493
1494	ncell=0
1495
1496	IF (Field(ind_b)%field_type==field_T) THEN
1497	nvert=6
1498
1499	DO ind=ind_b,ind_e
1500	d=>domain(ind)
1501
1502	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1503	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1504	IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
1505	ENDDO
1506	ENDDO
1507
1508	END DO
1509
1510	CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
1511
1512	displ=0
1513	DO i=1,mpi_rank
1514	displ=displ+ncell_glo(i-1)
1515	ENDDO
1516	FieldVarId(NbField)%displ=displ
1517	ncell_tot=sum(ncell_glo(:))
1518
1519	! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
1520	FieldId(NbField)=ncid
1521	status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
1522	status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1523
1524	IF (Field(ind_b)%ndim==2) THEN
1525	FieldVarId(NbField)%size=1
1526	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1527	ELSE IF (Field(ind_b)%ndim==3) THEN
1528	FieldVarId(NbField)%size=1
1529	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1530	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1531	ELSE IF (Field(1)%ndim==4) THEN
1532	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1533	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1534	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1535	! status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
1536	ENDIF
1537
1538	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1539
1540	status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1541	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1542	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1543	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1544	status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1545	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1546	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1547	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1548	status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1549	status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1550
1551	IF (Field(ind_b)%ndim==2) THEN
1552	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1553	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1554	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
1555	ELSE IF (Field(ind_b)%ndim==3) THEN
1556	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1557	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1558	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, &
1559	(/ncell_tot,size(field(ind_b)%rval3d,2),1/))
1560	ELSE IF (Field(ind_b)%ndim==4) THEN
1561	DO i=1,FieldVarId(NbField)%size
1562	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /), &
1563	FieldVarId(NbField)%nc_id(i))
1564	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
1565	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
1566	(/ncell_tot,size(field(ind_b)%rval4d,2),1/))
1567	ENDDO
1568	ENDIF
1569
1570
1571	status = NF90_ENDDEF(ncid)
1572
1573	ncell=1
1574	DO ind=ind_b,ind_e
1575	d=>domain(ind)
1576
1577	n=0
1578	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1579	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1580	IF (single .OR. d%own(i,j)) n=n+1
1581	ENDDO
1582	ENDDO
1583
1584	ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1585
1586	n=0
1587	DO j=d%jj_begin-halo_size,d%jj_end+halo_size
1588	DO i=d%ii_begin-halo_size,d%ii_end+halo_size
1589	IF (d%own(i,j) .OR. single) THEN
1590	n=n+1
1591	CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
1592	lon(n)=lon(n)*180/Pi
1593	lat(n)=lat(n)*180/Pi
1594	DO k=0,5
1595	CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
1596	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1597	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1598	ENDDO
1599	ENDIF
1600	ENDDO
1601	ENDDO
1602	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
1603	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
1604	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1605	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1606
1607	ncell=ncell+n
1608	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1609	END DO
1610
1611	ELSE IF (Field(ind_b)%field_type==field_Z) THEN
1612	nvert=3
1613	DO ind=ind_b,ind_e
1614	d=>domain(ind)
1615
1616	DO j=d%jj_begin+1,d%jj_end
1617	DO i=d%ii_begin,d%ii_end-1
1618	ncell=ncell+1
1619	ENDDO
1620	ENDDO
1621
1622	DO j=d%jj_begin,d%jj_end-1
1623	DO i=d%ii_begin+1,d%ii_end
1624	ncell=ncell+1
1625	ENDDO
1626	ENDDO
1627
1628	END DO
1629
1630	CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
1631
1632	displ=0
1633	DO i=1,mpi_rank
1634	displ=displ+ncell_glo(i-1)
1635	ENDDO
1636	FieldVarId(NbField)%displ=displ
1637	ncell_tot=sum(ncell_glo(:))
1638
1639	! status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
1640	! status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
1641	FieldId(NbField)=ncid
1642	status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
1643	status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
1644
1645	IF (Field(ind_b)%ndim==2) THEN
1646	FieldVarId(NbField)%size=1
1647	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1648	ELSE IF (Field(ind_b)%ndim==3) THEN
1649	FieldVarId(NbField)%size=1
1650	ALLOCATE(FieldVarId(NbField)%nc_id(1))
1651	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
1652	ELSE IF (Field(1)%ndim==4) THEN
1653	FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
1654	ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
1655	status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
1656	ENDIF
1657
1658
1659
1660	status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
1661
1662	status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
1663	status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
1664	status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
1665	status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
1666	status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
1667	status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
1668	status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
1669	status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
1670	status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
1671	status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
1672
1673
1674	IF (Field(ind_b)%ndim==2) THEN
1675	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
1676	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1677	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
1678	ELSE IF (Field(ind_b)%ndim==3) THEN
1679	status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
1680	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
1681	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, &
1682	(/ncell_tot,size(field(ind_b)%rval3d,2),1/))
1683	ELSE IF (Field(ind_b)%ndim==4) THEN
1684	DO q=1,FieldVarId(NbField)%size
1685	status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec, &
1686	(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
1687	status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
1688	status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
1689	(/ncell_tot,size(field(ind_b)%rval4d,2),1/))
1690	ENDDO
1691	ENDIF
1692
1693	status = NF90_ENDDEF(ncid)
1694
1695	ncell=1
1696	DO ind=ind_b,ind_e
1697	d=>domain(ind)
1698	CALL swap_geometry(ind)
1699	CALL swap_dimensions(ind)
1700
1701	n=0
1702	DO j=jj_begin+1,jj_end
1703	DO i=ii_begin,ii_end-1
1704	n=n+1
1705	ENDDO
1706	ENDDO
1707
1708	DO j=jj_begin,jj_end-1
1709	DO i=ii_begin+1,ii_end
1710	n=n+1
1711	ENDDO
1712	ENDDO
1713
1714	ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
1715
1716	n=0
1717
1718	DO j=jj_begin+1,jj_end
1719	DO i=ii_begin,ii_end-1
1720	nij=(j-1)*iim+i
1721	n=n+1
1722	CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
1723	lon(n)=lon(n)*180/Pi
1724	lat(n)=lat(n)*180/Pi
1725	CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1726	CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1727	CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1728
1729	DO k=0,2
1730	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1731	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1732	ENDDO
1733	ENDDO
1734	ENDDO
1735
1736	DO j=jj_begin,jj_end-1
1737	DO i=ii_begin+1,ii_end
1738	nij=(j-1)*iim+i
1739	n=n+1
1740	CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
1741	lon(n)=lon(n)*180/Pi
1742	lat(n)=lat(n)*180/Pi
1743	CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
1744	CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
1745	CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
1746
1747	DO k=0,2
1748	bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
1749	bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
1750	ENDDO
1751	ENDDO
1752	ENDDO
1753
1754
1755	status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
1756	status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
1757	status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1758	status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
1759
1760	ncell=ncell+n
1761	DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
1762	END DO
1763	ENDIF
1764
1765
1766	END SUBROUTINE Create_Header_mpi
1767
1768	SUBROUTINE Close_files
1769	USE netcdf
1770	IMPLICIT NONE
1771	INTEGER :: i,k,status
1772	!$OMP MASTER
1773	DO i=1,NbField
1774	status=NF90_CLOSE(FieldId(i))
1775	ENDDO
1776	!$OMP END MASTER
1777	END SUBROUTINE Close_files
1778
1779
1780	function int2str(int)
1781	implicit none
1782	integer, parameter :: MaxLen=10
1783	integer,intent(in) :: int
1784	character(len=MaxLen) :: int2str
1785	logical :: flag
1786	integer :: i
1787	flag=.true.
1788
1789	i=int
1790
1791	int2str=''
1792	do while (flag)
1793	int2str=CHAR(MOD(i,10)+48)//int2str
1794	i=i/10
1795	if (i==0) flag=.false.
1796	enddo
1797	end function int2str
1798
1799	end module write_field
1800

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

Download in other formats: