Context Navigation

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

Last change on this file since 82 was 82, checked in by ymipsl, 12 years ago

improvment of the output and CF conformity

YM

File size: 62.1 KB

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

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

Download in other formats: