Context Navigation

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

Last change on this file since 378 was 378, checked in by ymipsl, 8 years ago

Add vertical axis definition in output field (merged from heat/aquaplanet).

YM

File size: 63.3 KB

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

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

Download in other formats: